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conditional density of given X n _^ is derived. This 
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variate extensions of the model are considered briefly. 

Second, three methods for generating first-order 
moving average sequences with Exponential marginals are 
examined. These generalize the EMA(l) Exponential model. 
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gated in the moving average models. 

A preliminary analysis of wind speed data obtained over 
a 15 year period in the Gulf of Alaska is presented. A 
model with four harmonic deterministic mean multiplying 
random innovative factors modeled by a GLAR(l) process is 
developed. Correlograms and periodograms are used to deter- 
mine the model for the mean and the structure of the 
innovation process. 



4 



TABLE OF CONTENTS 



I. INTRODUCTION 

II. MODELS WITH GAMMA MARGINALS 

A. INTRODUCTION 

B. FIRST-ORDER AUTOREGRESSIVE BETA-GAMMA 

MODEL, GLAR(l) 

1. Introduction 

2. Correlation Structure 

3. Conditional Expectation and Conditional 

Density 

4 . Maximum Likelihood Estimation 

5 . Directional Moments 



C. FIRST-ORDER MOVING AVERAGE BETA GAMMA MODEL - 

1 . Introduction 

2. Correlation Structure 

3. Directional Moments 

4 . Empirical P (X n+ ^ > X ) 

D. OTHER CASES OF THE GLARMA ( P , Q ) MODEL 

1 . Introduction 

2. GLARMA (1,1) 

3 . GLAR (2) 

4. BGAR(l), BIVARIATE MODEL 

E. NUMERICAL CONVERGENCE OF THE MAXIMUM LIKE- 

LIHOOD COMPUTER PROGRAM AND SIMULATION 
STUDY OF PROPERTIES OF ESTIMATORS 

III. MOVING AVERAGE MODELS 

A. INTRODUCTION 



19 

23 

23 

26 

26 

29 

32 

36 

45 

47 

50 

50 

51 

60 

62 

62 

62 

63 

63 

77 

35 

99 

99 



5 



1 . Negative Correlation 



101 



B. 

C. 



D. 



E. 



2. New Exponential Moving Average Model 

of Order One, NGMA(l) 

3 . The Moving Minimum Model 

4 . 'The Beta-Exponential Model 



NEGATIVE CORRELATION IN MOVING AVERAGE 
MODELS 

THE EXTENDED EMA(l) MODEL, NEMA(l) 



1 . 

2 . 

3. 

4 . 
5. 

6 . 
7. 

THE 

1 . 

2 . 

3. 

4 . 

5. 

6 . 



Introduction 

Correlation Structure 



P (X 



n+1 



V 



Laplace Transform of Sums 



Laplace Transform of the Distribution 
of Counts 



The Spectrum of Counts 

Joint Laplace-Stielt j es Transform 

of X and X ,, 

n n+1 

MOVING MINIMUM MODEL 



Introduction 

Correlation Structure 

Negative Correlation 

Joint Density of X p and X n+ ^ 

Conditional Expectation and P(X^ + -^>X n ) “ 

Conditional Expectation and P(X , > X ) 
for Antithetic Variables 



THE BETA-EXPONENTIAL MODEL 

1. Introduction 

2. Correlation Structure, Positive and 

Negative 



102 

102 

103 

103 

108 

108 

111 

120 

127 

143 

147 

156 

159 

159 

161 

164 

168 

179 

193 

209 

209 

210 



6 



3. Laplace-Stielt j es Transform of a Sum 213 

4. Empirical P (X n+ ^ > X ) 214 

IV. PRELIMINARY DATA ANALYSIS 228 

A. INTRODUCTION 22 8 

B. ANALYSIS OF THE RAW DATA 229 

C. THE FORM OF THE MEAN; DETRENDING THE DATA 255 

D. RESIDUAL PROCESS PROBABILITY STRUCTURE 282 

E. REFINING THE FORM OF THE MEAN; A FURTHER 

DETRENDING 307 

F. RESIDUAL ANALYSIS 331 

G. A FURTHER REFINEMENT OF THE MEAN; THE FINAL 

DETRENDING 339 

H. SUMMARY 347 

LIST OF REFERENCES 352 

INITIAL DISTRIBUTION LIST 355 



7 



LIST OF TABLES 



II.E.l Parameter Values for the Simulation Study - 86 

II. E. 2 Results of Search for Maximum Likelihood 

Estimates in a GLAR(l) Model; k = 4.0, 
q = 1.0 88 

II. E. 3 Results of Search for Maximum Likelihood 

Estimates in a GLAR(l) Model; k = 4.0, 
q = 3.0 92 

II. E. 4 Results of Search for Maximum Likelihood 

Estimates in a GLAR(l) Model; k = 0.75, 
q = 0.1875 95 

II. E. 5 Summary of Simulation Results 98 

III. C. 3.1 Formulae for P(X . > X ) for the 

NEMA(l) Model — 125 

III. C. 3. 2 Computations of P (X . > X ) for the 

NEMA(l) Model 12 6 

IV. B.l Average Wind Velocity by Month 248 

IV. B. 2 Estimated Correlations for Raw Data 253 

IV.C.l Estimated Coefficients for Models of 

the Mean 257 

IV.D.l Moment Estimate of the Gamma Shape 

Parameter by Year 29 9 

IV. D. 2 Effect of Cuts on Moment Estimates 300 

IV. D. 3 Estimated Correlations for 1 Harmonic 

Detrended Data 302 

IV.E.l Estimated Correlations for 2 Harmonic 

Detrended Data 332 

IV.G.l Estimated Correlations for 4 Harmonic 

Detrended Data 349 



8 



LIST OF FIGURES 



II. B. 6.1 GLAR(l) Sample Path; k = 4.0; q = 1.0 48 

II.B.6.2 GLAR(l) Sample Path; k = 4.0; q = 3.0 49 

II. D. 2.1 GLARMA (1,1) Model; Rho(l) vs Rho(2) 69 

II. D. 3.1 GLAR(2) Model; Rho(l) vs Rho(2) 76 

II.E.l M.L.E. and Moment Estimates; 

TRUE k = 4.0; TRUE q = 1.0 90 

II. E. 2 M.L.E. and Moment Estimates; 

TRUE k = 4.0; TRUE q = 3.0 93 

II. E. 3 M.L.E. and Moment Estimates; 

TRUE k = 0.75; TRUE q = 0.18 75 96 

II I. C. 3.1 NEMA Sample Path; Alpha: 0.95 

Beta: 0.50; TRUE Rho : 0.25 12 8 

III. C. 3. 2 NEMA Sample Path; Alpha: 0.70 

Beta: 0.95; TRUE Rho: 0.22 129 

III. C. 3.3 NEMA Sample Path; Alpha: 0.30 

Beta: 0.90; TRUE Rho: 0.20 130 

III. C. 3. 4 NEMA Scatter Plot; Alpha: 0.95 

Beta: 0.50; TRUE Rho: 0.25 131 

III. C. 3. 5 NEMA Scatter Plot; Alpha: 0.70 

Beta: 0.95; TRUE Rho: 0.22 132 

III. C. 3. 6 NEMA Scatter Plot; Alpha: 0.30 

Beta: 0.90; TRUE Rho: 0.20 13 3 

II I. C. 3. 7 NEMA Sample Path; Alpha: 0.95 

Beta: 0.50; TRUE Rho : -0.16 134 

III. C. 3. 8 NEMA Sample Path; Alpha: 0.70 

Beta: 0.95; TRUE Rho: -0.14 135 

III. C. 3. 9 NEMA Sample Path; Alpha: 0.30 

Beta: 0.90; TRUE Rho: -0.13 136 

III. C. 3. 10 NEMA Scatter Plot; Alpha: 0.95 

Beta: 0.50; TRUE Rho: -0.16 13 7 



9 



III. C. 3. 11 NEMA Scatter Plot; Alpha: 0.70 

Beta: 0.95; TRUE Rho : -0.14 138 

III. C. 3. 12 NEMA Scatter Plot; Alpha: 0.30 

Beta: 0.90; TRUE Rho: -0.13 139 

III. C. 6.1 NEMA Spectrum of Counts 153 

III. C. 6. 2 NEMA Spectrum of Counts - 

Constant P(X , . > X ) 154 

n+1 n 

III. C. 6. 3 NEMA Spectrum of Counts - 

Constant Correlation 155 

III. D. 3.1 Moving Minimum Model - Range of 

Values for Positive Correlation 169 

II I. D. 3. 2 Moving Minimum Model - Range of 

Values for Negative Correlation 170 

III. D. 4.1 Regions for Values of Z Given Y > bX 172 

III. D. 4. 2 Regions for Values of Z Given Y £ bX 172 

II I. D. 5.1 Moving Minimum Model - P(X . > X ): 

Positive Correlation 185 

III. D. 5. 2 Min Model Sample Path - B Value: 1.00 — 187 

III. D. 5. 3 Min Model Sample Path - B Value: 1.50 — 188 

III. D. 5. 4 Min Model Sample Path - B Value: 0.67 — 189 

III. D. 5. 5 Min Scatter Plot - B Value: 1.00 190 

III. D. 5. 6 Min Scatter Plot - B Value: 1.50 191 

III. D. 5. 7 Min Scatter Plot - B Value: 0.67 192 

III. D. 6.1 Moving Minimum Model - P (X . > X ): 

Negative Correlation 20 2 

III. D. 6. 2 Min Model Sample Path - B Value: 1.00 — 203 

III. D. 6. 3 Min Model Sample Path - B Value: 1.50 — 204 

III. D. 6. 4 Min Model Sample Path - B Value: 0.67 — 205 

III. D. 6. 5 Min Scatter Plot - B Value: 1.00 206 

III.D.6.6 Min Scatter Plot - B Value: 1.50 207 



10 



III. D. 6. 7 Min Scatter Plot - B Value: 0.67 208 

III. E. 4.1 Beta Sample Path - q: 0.25 - 

True Rho: 0.19 216 

III. E. 4. 2 Beta Sample Path - q: 0.67 - 

True Rho: 0.22 217 

III. E. 4. 3 Beta Sample Path - q: 0.33 - 

True Rho: 0.22 218 

III. E. 4. 4 Beta Scatter Plot - q: 0.25 - 

True Rho: 0.19 219 

III. E. 4. 5 Beta Scatter Plot - q: 0.67 - 

True Rho: 0.22 220 

III. E. 4. 6 Beta Scatter Plot - q: 0.33 - 

True Rho: 0.22 221 

III. E. 4. 7 Beta Sample Path - q: 0.25 - 

True Rho: -0.12 222 

III. E. 4. 8 Beta Sample Path - q: 0.67 - 

True Rho: -0.14 223 

III. E. 4. 9 Beta Sample Path - q: 0.33 - 

True Rho: -0.14 224 

III. E. 4. 10 Beta Scatter Plot - q: 0.25 - 

True Rho: -0.12 225 

III. E. 4. 11 Beta Scatter Plot - q: 0.67 - 

True Rho: -0.14 226 

III. E. 4. 12 Beta Scatter Plot - q: 0.33 - 

True Rho: -0.14 22 7 

.. IV.B.la 19 55 Raw data 230 

IV. 3. lb 19 56 Raw data 231 

IV.B.lc 1957 Raw data 232 

IV.B.ld 1958 Raw data 233 

IV.B.le 1959 Raw data 234 

IV.B.lf 1960 Raw data 235 

IV.B.lg 1961 Raw data 236 



11 



237 



IV.B.lh 1962 Raw data 

IV.B.li 19 6 3 Raw data 238 

IV.B.lj 1964 Raw data 239 

IV.B.lk 19 6 5 Raw data 24 0 

IV.B.li 1966 Raw data 241 

IV.B.lm 19 6 7 Raw data 24 2 

IV.B.ln 19 6 8 Raw data 243 

IV.B.lo 1969 Raw data 24 4 

IV.B.lp 15 year average raw data 245 

IV. B. 2 Monthly average wind speed 24 7 

IV. B. 3 Periodogram, 15 year average raw data 250 

IV. B. 4 Log periodogram, 15 year average 

raw data 251 

IV.C.l Average data plotted against 

1 harmonic smoothed data 259 

IV. C . 2a 19 55 Detrended data 261 

IV.C.2b 1956 Detrended data 262 

IV.C.2c 1957 Detrended data 263 

IV.C.2d 1958 Detrended data 264 

IV.C.2e 1959 Detrended data 265 

.. IV.C.2f 1960 Detrended data 266 

IV.C.2g 1961 Detrended data 267 

IV.C.2h 1962 Detrended data 268 

IV. C . 2i 1963 Detrended data 269 

IV.C.2j 1964 Detrended data 2/0 

IV.C.2k 1965 Detrended data 271 

IV.C.21 1966 Detrended data 272 



12 



273 



IV.C.2m 1967 Detrended data 

IV.C.2n 1968 Detrended data 274 

IV.C.2o 19 69 Detrended data 275 

lV.C.2p 15 year average detrended data 

(exponential sine) 276 

IV. C. 3 15 year average detrended data (sine) — 277 

IV. C. 4 Log periodogram, 15 year average 

detrended (sine) data 278 

IV. C. 5 Five step moving average log 

periodogram, 15 year average 

detrended (sine) data 279 

IV. C. 6 Log periodogram, 15 year average 

detrended (exponential sine) data 280 

IV. C. 7 Five step moving average log 

periodogram, 15 year average 

detrended (exponential sine) data 281 

IV.D.la Histogram and box plot, 

19 55 detrended data 28 3 

IV.D.lb Histogram and box plot, 

1956 detrended data 284 

IV.D.lc Histogram and box plot, 

19 57 detrended data 285 

IV.D.ld Histogram and box plot, 

1958 detrended data 286 

IV.D.le Histogram and box plot, 

1959 detrended data 287 

IV.D.lf Histogram and box plot, 

1960 detrended data 288 

IV.D.lg Histogram and box plot, 

19 61 detrended data 2 89 

IV.D.lh Histogram and box plot, 

1962 detrended data 29 0 

IV.D.li Histogram and box plot, 

196 3 detrended data 291 



13 



IV.D.lj Histogram and box plot, 

1964 detrended data 292 

IV.D.lk Histogram and box plot, 

1965 detrended data 293 

IV.D.ll Histogram and box plot, 

1966 detrended data 294 

IV.D.lm Histogram and box plot, 

1967 detrended data 295 

IV.D.ln Histogram and box plot, 

196 8 detrended data 296 

IV.D.lo Histogram and box plot, 

1969 detrended data 297 

IV.D.lp Histogram and box plot, pooled 

data from 1955-1969 298 

IV. D. 2 Periodogram 1955 detrended data 

with AR (1) spectrum 303 

IV. D. 3 Log periodogram 1955 detrended data 

with log AR ( 1 ) spectrum 304 

IV. D. 4 Periodogram 1969 detrended data 

with AR(1) spectrum 305 

IV. D. 5 Log periodogram 1969 detrended data 

with log AR(1) spectrum 306 

IV.D.6a Correlogram 1955 detrended data 308 

IV.D.6b Correlogram 1956 detrended data 309 

IV.D.6c Correlogram 1957 detrended data 310 

IV.D.6d Correlogram 1958 detrended data 311 

IV.D.6e Correlogram 1959 detrended data 312 

IV.D.6f Correlogram 1960 detrended data 313 

IV.D.6g Correlogram 1961 detrended data 314 

IV.D.6h Correlogram 1962 detrended data 315 

IV.D.6i Correlogram 1963 detrended data 316 

IV. D . 6 j Correlogram 1964 detrended data 317 



14 



IV.D. 6 k Correlogram 1965 detrended data 318 

IV.D.61 Correlogram 1966 detrended data 319 

IV.D. 6 m Correlogram 1967 detrended data 320 

IV.D. 6 n Correlogram 1968 detrended data 321 

IV.D .60 Correlogram 1969 detrended data 322 

IV.D. 6 p Average correlogram, detrended data 323 

IV. D. 7 Correlogram of average detrended data 324 

IV.E.l Average data plotted against 2 

harmonic smoothed data 326 

IV. E. 2 Periodogram 1955 2 harmonic detrended 

data with AR(1) spectrum 327 

IV. E. 3 Log periodogram 1955 2 harmonic 

detrended data with log AR(1) spectrum — 328 

IV. E. 4 Periodogram 1969 2 harmonic detrended 

data with AR(1) spectrum 329 

IV. E. 5 Log periodogram 1969 2 harmonic 

detrended data with log AR(1) spectrum — 330 

IV. E . 6 Periodogram 15 year average 2 

harmonic detrended data 333 

IV. E. 7 Log periodogram 15 year average 

2 harmonic detrended data 334 

IV.F.l Periodogram prewhitened 15 year 

average 2 harmonic detrended data 336 

IV. F. 2 Log periodogram prewhitened 15 year 

average 2 harmonic detrended data 337 

IV.F.3 Correlogram prewhitened 15 year average 

2 harmonic detrended data 338 

IV.G.l Periodogram 1955 4 harmonic detrended 

data with AR(1) spectrum 341 

IV. G. 2 Log periodogram 1955 4 harmonic 

detrended data with log AR(1) spectrum — 342 



15 



IV. G. 3 Periodogram 1969 4 harmonic detrended 

data with AR(1) spectrum 343 

IV. G. 4 Log periodogram 1969 4 harmonic 

detrended data with log AR(1) spectrum — 344 

IV. G. 5 Periodogram 15 year average 4 harmonic 

detrended data 345 

IV. G. 6 Log periodogram 15 year average 

4 harmonic detrended data 346 

IV. G. 7 GLAR(l) Sample path; k = 4.0, q = 0.6; 

TRUE Rho: 0.85 34 8 



16 



ACKNOWLEDGEMENT 



I want to take this opportunity to express my gratitude 
to the following individuals: 

Professor P. A. W. Lewis whose efficiency in getting 
my research started and whose unflagging enthusiasm and 
support are largely responsible for getting this project 
completed . 

Professor R. R. Read whose counsel and personal support 
were as valuable as his technical knowledge when I was 
preparing for my qualifying exams. 

Associate Professor J. K. Hartman who helped even when 
he couldn't afford the time and whose suggestion concerning 
the search technique for the maximum likelihood program 
saved many hours of CPU time. 

Associate Professor F. R. Richards who filled the gap 
with virtually no notice. 

Professors G. H. Bradley and R. J. Renard who took 
time from busy schedules to help an outsider. 

Professor D. R. Barr whose wise advice may have made 
the difference between success and failure. 

Distinguished Professor Emiritus F. D. Faulkner and 
Associate Professor R. H. Franke who showed me how to 
eliminate the singularities in a conditional density. 

Assistant Professor R. W. Garwood who provided the 
data for analysis and lent his oceanographical expertise 
to the evaluation of the analysis. 



17 



Ladies and Gentlemen of W. R. Church Computer Center, 
especially those on the 1600-2400 shift, whose support 
was invaluable during periods when plotter and CPU use 
were intense - 

Nancy Ritz Hugus whose support and understanding make 
all good things possible. 



18 



INTRODUCTION 



I . 

The classical approach to time series analysis based on 
linear, additive models with normally distributed, constant 
variance residuals is probably best presented by the work of 
Box and Jenkins [Ref. 1] . Although their work is widely ac- 
cepted and used, it is not applicable to some important time 
series. This is mainly because the Box-Jenkins approach is 
based on an assumed normal distribution for the series in 
question. However, the assumption of normality is not appro- 
priate when the series is known to be non-negative. Such 
series typically involve times between successive events in 
event processes. Examples are easy to construct. Times be- 
tween arrivals at a hospital emergency room, times between 
breakdowns in a tank main drive assembly, and times between 
detections of enemy armor vehicles are a sample of series of 
this type. Because of the non-negative nature of the series, 
the Box-Jenkins distributional assumptions and, hence, the 
analysis techniques are inappropriate. There is, of course, 
the possibility of data transformations but this is not appro- 
priate with very skewed marginal distributions and it is, in 
most cases, difficult to ascertain what the transformation 
does to the correlation structure of the series. 

Gaver and Lewis [Ref. 2] wrote the pioneering paper on 
the subject of autoregressive processes with non-normal marginal 
distributions. They presented the method for determining the 
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distribution of the innovative terms in the basic, linear, 
additive, autoregressive equations (first-order stochastic 
different equation) 



X 



n 



PX 



n-1 



+ 



e 

n 



( 1 . 1 ) 



that was required to produce a given marginal distribution 

for the (X } sequence. They presented results for (X } se- 
n n 

quences with Exponential, Gamma, and mixed Exponential 
marginals. They also showed that this problem was the same 
as that of determining the class of self-decomposable (Type 
L) random variables (Feller, [Ref. 3], Loeve [Ref. 4]) al- 
though the connection between the solution to the self-decom- 
posable problem and equation (1.1) was not explicit in the 
literature. 

The Gaver and Lewis paper was followed by other papers 
which extended these results. Lawrance and Lewis [Ref. 5] 
presented a first-order moving average process with Exponential 
marginals. Jacobs and Lewis [Ref. 6] propounded a mixed auto- 
regressive-moving average of order one, EARMA(1,1), and 
Exponential marginals. This was extended to an arbitrary 
order EARMA(p,q) process by Lawrance and Lewis [Ref. 7]. A 
further refinement of the first-order. Exponential, auto- 
regressive process (NEAR(l)) was presented by Lawrance and 
Lewis [Ref. 8] . While this contained the previous EAR ( 1 ) 
model, it did not suffer from the degeneracy inherent in (1.1) . 
Jacbos applied these models to closed cyclic queueing networks 
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[Ref. 9] and Lewis and Shedler applied them to models of 
computer processes [Ref. 10] . 

This paper extends the results of these researchers and 
others in three areas. In Chapter II a mixed ARMA(p,q) model 
with Gamma marginals proposed by Lewis [Ref. 11] , the GLARMA(p,q) 
model, is examined. The correlation structure is derived for 
several values of p and q. Of particular note is the AR(1) 
case (p = 1, q = 0) , called GLAR(l) , where the conditional 
density of X n given X n _^ is derived. This leads to the deri- 
vation of a likelihood function and a numerical technique to 
evaluate and maximize the likelihood function with respect to 
the model parameters. This provides a useful technique for 
estimating model parameters. Using this numerical technique 
a simulation study of the properties of maximum likelihood 
estimators for the parameters of the model is given. 

The correlation structure is derived for other models in 
the GLARMA(p,q) family: the first-order moving average, the 

second-order autoregressive, the first-order mixed autoregres- 
sive-moving average and a bivariate first-order autoregressive 
process. These different models, particularly the bivariate 
extension, demonstrate the flexibility of the GLARMA(p,q) 
model . 

In Chapter III the first-order moving average process with 
Exponential marginals of Lawrance and Lewis [Ref. 5] is ex- 
tended to a two parameter model. This is done by utilizing 
the NEAR ( 1 ) structure which combines two independent Exponen- 
tial random variables into a random variable with Exponential 
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distribution. A fairly complete set of characteristics of 
this model are derived. In particular the correlation struc- 
ture, the quantity P (X n+ ^ > X ) , the Laplace transform of 
sums of X n 's, the Laplace transforms of the distribution of 
counts, the (Bartlett) spectrum of counts, and the joint 
Laplace-Stielt jes transforms of X n and X n+1 are addressed. 

These characteristics are compared to those of other proc- 
esses which produce marginally Exponential random variables. 

In Chapter IV the models of Chapter II are used in a 
preliminary data analysis of wind speed data. This repre- 
sents the first effort to apply these models to a large, real 
world data base. A model for simulating wind data is pre- 
sented and parameter estimates for the data are derived using 
the numerical approximation to the maximum likelihood presented 
in Chapter II. 
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II. MODELS WITH GAMMA MARGINALS 



A. INTRODUCTION 

There have been several schemes suggested for the modeling 
of dependent random variables with Gamma marginals. The 
Gamma autoregressive process of order one (GAR(l) ) by Gaver 
and Lewis [Ref. 2], the discrete autoregressive process of 
order one (GDAR(l)) by Jacobs and Lewis [Ref. 12 ], the Gamma 
Beta autoregressive process of order one (GBAR(l)) by Fishman 
[Ref.. 13] and Lawrance and Lewis [Ref. 14], and the Gamma auto- 
regressive process of order one (GLAR(i)) by Lewis [Ref. 10]. 
There is also an attempt to use multivariate Gammas obtained 
by the inverse probability integral transform in a time series 
context by Schmeiser [Ref. 15] . 

The GAR(l) model generates an {X^} series using the stan- 
dard first-order autoregression equation (first-order stochas- 
tic difference equation) 

X = pX . + e , 0 < p < 1. (II. A. 1) 

n H n-1 n — 

The innovative factor, e , has Laplace-Stielt jes transform 

\ k. 

of (p + (1— p) "t~ — ) and the {X } random variables have marginal 

A ' s n. 

Gamma distributions with shape parameter k and scale parameter 
A. The marginal density function of the (X n ) random variable 
is 
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f x (x;k,A) 



k-1 -Ax 
x e 



(II .A. 2) 



f x ( x ; k , A ) 
n 



x' ' ' ' rck) 

A > 0, x > 0, k > 0. 



The model tends to produce runs of decreasing value when inno- 
vative term has successive realizations of value zero. The 
GAR(l) is in this sense highly degenerate, even though it is 
a true linear process. Ad hoc estimates of model parameters 
are available which produce the exact p value if the series 
is long enough [Ref. 2] . However, maximum likelihood esti- 
mates have not been produced. This model is not extendable 
to a moving average process. 

The GDAR(l) produces an {X n } sequence using the first- 
order autoregressive equation with random coefficients. 



X 



n 



X . 

n n-1 



(1 " V n )G 



n' 



(II .A. 3) 



where {V n , n = 1,2,...} is an iid sequence of binary random 

variables with P(V = 1) = 1-P(V =0) = p, (G , n = 1,2,...} 

n n n 

is an iid Gamma sequence. 

This sequence produces runs of constant value when suc- 
cessive realizations for V produce value 1. When V equals 

n n 

zero, a new value is selected. Obviously, this model has 
limited value in general applications and is even more degener- 
ate than GAR(l) process. 

The GBAR (1) is the most flexible model in that it contains 
the GAR ( 1) and GLAR(l) models as special cases. It produces 
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-? 



an (X n > sequence using 



X = BB X 
n n n+1 



n' 



(II .A. 4) 



where (B n , n = 1,2,...} is an iid Beta (k-q,q) sequence. e n 
was shown by Lawrance and Lewis to be the sum of a Gamma 
variable and the innovation process of the GAR(l) process 
[Ref. 14] . Although flexible in the sense that it contains the 
other models, it can not be extended to a moving average proc- 
ess. In addition, conditional densities and, hence, maximum 
likelihood estimates are not available. This is because the 
innovation random variable for the GAR(l) process, while it 
can be generated as a random sum of random variables, does not 
have a known distribution function. 

The most valuable and flexible model seems to be the 
GLAR(l) which produces an {X n } sequence using the stochastic 
difference equation with random coefficients 



X 



n 



= B X . 
n n-1 



C G , 
n n' 



(II .A. 5) 



where {X , n = 0,1,...} is a second-order stationary sequence 
n 

of Gamma random variables, (B^, n = 1/2,...} and {C n , n = 1,2,..._ 
are iid Beta random variables, and {G n } is an iid sequence of 
Gamma random variables. This model is extendable as an auto- 
regressive process of arbitrary order and as a moving average 
process of arbitrary order. These two forms can also be combined 
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to form a mixed model, the so-called Gamma Lewis autoregressive - 
moving average process of orders p and q (GLARMA (p, q) ) . 

This chapter of the thesis examines some of the special 
cases of the model. One case in particular, the AR(1) form, 
is reasonably extensively examined. The correlation struc- 
ture is developed. The conditional expectation and density 
are derived. The latter is used as the basis of a numerical 
approximation to the maximum likelihood method of parameter 
estimation. Directional moments and the probability of X n+ ^ 
being greater than X n are derived. In a later cnapter of 
this thesis, this model is used as a basis for analyzing wind 
speed data. 

The special case of the moving average of order one is 
examined in seme detail. The correlation structure is derived 
with some emphasis on exploring the restrictions on the range 
of correlations that are possible. Directional moments and 
an empirical examination of the probability that X n+ ^ is 
greater than X n are examined. 

As a demonstration of the flexibility and extendability 
of the model the mixed model of order one, the autoregressive 
model of order two, and a bivariate model are introduced and 
their correlation structures derived. 

B. FIRST-ORDER AUTOREGRESSIVE BETA-GAMMA MODEL, GLAR(l) 

1 . Introduction 

The first-order autoregressive Beta-Gamma model is a 
special case of the GLARMA(p,q) model when q = 0 and p = 1. 
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The autoregressive model generates an {X n } sequence using 



X 



n 



= AX , + 
n n-1 



B G , 
n n 



(II. B. 1.1) 



where (X n , n = 0,1,...} is a second-order stationary sequence 

of random variables with a Gamma (k,l) marginal distribution; 

{A n , n = 1,2,...} is an iid sequence of Beta(k-q,q) random 

variables; {B n , n = 1,2,...} is an iid sequence of Beta(q,k-q) 

random variables; (G n , n = 0,1,...} is an iid sequence of 

Gamma (k,l) random variables; {A } , {B } , and (G } are inde- 

n n n 

pendent; 0 < q £ k. 

Choosing Xq = Gq makes the {X n } sequence stationary. 

The Gamma random variables are parameterized by the shape 
parameter and the mean, rather than the scale parameter. This 
somewhat unusual parameterization has some advantages in 
statistical work since Gamma(k,y) = y Gamma(k,l) [Ref. 7], 

A Gamma (k, 1) random variable has density 

f G ( x ; k , 1 ) = YTkT xk_1 e " kX ' x > °' k > 0 (II. B. 1.2) 



This is a special case of the more general density 

,k. k kx 

f G (x.-k,u) = x k '! e 11 

where y = 1. Of course, since the scale parameter, shape 

]c 

parameter and mean are related by the relation y = V-, the den- 
sity can be specified by any two of the three parameters. The 
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typical parameterization in terms of the scale parameter, X, 
is useful because G(k 1 ,A) + G(k 2 ,A) = G(k 1 +k 2 ,X). This re- 
lationship is not true when the parameterization is through 
the shape parameter, k, and the mean, y. 

A Beta(m,n) random variable has density 



f 0 (x;m,n) 



0 



r (m+n) 

r(m) T(n) 

< x < 



m-1 . n-1 

x (1-x) 

1 , m > 0 , 



/ 



n > 0 . 



(II. B. 1.3) 



For the Beta random variable to be properly defined 

each of the parameters must be positive. Hence, when q = k, 

(II. B. 1.1), as defined above, is no longer appropriate since 

each Beta random variable has a parameter that is identically 

zero. In this case when q = k it is understood that the (A } 

sequence is considered to be identically zero and the {B^} 

sequence one. Therefore, II. B. 1.1. becomes simply X n = G n , 

and the (X } seauence, like the (G } sequence, is iid. A 
n n 

justification of this generation scheme for a Gamma process 
as defined by II. B. 1.1 was provided by Lawrance and Lewis 
[Ref . 14 , pp . 24 ] . 

In this section the correlation structures of the J 'X n > 

sequence and that of the { X^ } and {G^} sequences are addressed. 

Other characteristics of the sequence, such as conditional 

expectation of X^ given X^_^, directional moments, and 

?(X , >X ) are considered. Of particular note is the dari- 

ng- 1 n 

vation of the conditional density of X n given X^_^. This 
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leads to the formulation of the likelihood function and a 
computer program to generate maximum likelihood estimates 
of parameter values. The numerical convergence properties 
of the likelihood method are assessed in Section II. E. 

2. Correlation Structure 



determined by a straightforward calculation. We have 



independent if j > i. Using these facts along with the iid 
nature and independence of the (B n ) and {G^} sequences yields 
the following expression when expectations are taken. 



The serial correlation of the {X n > series is easily 



X 



n 








E< V E(X n-l> + E(B „ )E(G n )E(x n-l 





Therefore 



C 0 V <V X n-l 




and 
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CORR (X , X .) 
n n- 1 



(II. B. 2.1) 



= = 1 - 0 < q < k. 

This is consistent with the fact that for q = k, {X r } is a 
sequence of i.i.d. Gamma variates which implies that 
CORR(X n X n _^) = 0. This correlation is easly extended by an 
induction argument to yield 

CORR (X ,X ) = (^) m , n > m > 0. (II.B.2.2) 

n n+m K — — 

The two sequences {X } and {G } can be viewed as a 

n n 

bivariate pair ( x n ' G n ) °f Gamma (k,l) random variables. 

Therefore, the correlation structure of these sequences may 
be of interest. Proceeding as before 



X 



n 



AX , + B G 
n n-1 n n 



X G 
n n 



AX .G + B G 
n n-1 n n n 



Taking expectations as before 



E (X G ) = E(A )E(X^ ,)E(G_) + E (B ) E (G ) 

n n n n-J. n n n 



e<x g ) = + 

n n k k ^2 



_ kza + k q + q 

~ k k 2 



E (X G ) = t+2 

n n k 2 



(II.B.2.3) 



Therefore 



COV (X , G ) = -% 

n n k 2 



and 
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(II .B.2 .4) 



CORR(X n ,G n > = | . 



When q = k this is 1 since X = G 

n n 



Pursuing the process one more step we determine 



CORR(X n ,G n _ 1 ) . 



X = AX . + B G , 
n n n-1 n n' 



XG , = AX .G ,+BGG , 

n n-1 n n-1 n-1 n n n-1 



Taking expectations as before and using II. B.2. 3 and the second- 
order stationarity of the (X n ) sequence, we get 



E(X„G 



n-1' 



Therefore, 



E (A ) E (X .C- ,) + E (B ) E (G ) E (G .) 

n n-1 n-i n n n-i 

( k-q ( k 2 +q , q.!.! = k 3 +kq-k 2 q-q 2 +qk 2 

1 k M ,2 ' + k .3 

k k 

k 3 +kq-q 2 



and 



COV(X ,G ,) = q(k 3 q) 

n n-1 k 3 

C0RE(X n' G n-l ) = ‘k 1 *^ 1 



(II .B.2 .5) 



(II .B.2. 6) 



II. B.2. 3 through II. B.2. 6 can be used in a simple induction 
argument to yield the general result 



CORR(X n ,G n _ m ) = (^ (^i^) 111 ' m= 0,1,..., n. (II. B.2. 7) 

When j is greater than i, X^^ and G ^ are independent. Hence, 
CORR(X i ,Gj) =0, j > i. 
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y is 



3 • Conditional Expectation and Conditional Density 
The conditional expectation of given X R ^ = 
trivially determined from the defining equation and the 
moments of the Beta distribution as 



E(X nl X n -l = y) " 




) y + 



a 



k* 



(II. B. 3.1) 



k cr 

Recognizing as the correlation between X and X , and 

x n n-1 

letting p be that correlation, II. B. 3.1 can be written as 



E(X J X n-l = y) = py + k = py + (1 ~ p) * (II.G.3.2) 

Thus the regression is linear in y. 

The conditional density of X n given X n _^ can also be 
determined. It is easiest to start by deriving the condi- 
tional distribution of X given X , = y. 

n 3 n-1 



plx n i x l K n -r') 



P (A y + B G < x) 
n n — 



P( 3 r. G ni x -V l 



Now [Ref. 14] the product of a Beta(q,k-q) random variable and 
a Gamma (k,l) random variable is a random variable with density 
Gamma (q, 3.) . Hence, if we let D n be a Gamma (q, 3) random 
variable. 



P(X < x I X , = y) = P (D < x-Av) (II.B.3.3) 

n — 1 n-1 n — n J 

This can be written as a convolution if one is careful about 
the upper limit of integration. Since D n is Gamma (q,^), it 
is non-negative. Hence, P(D n <_x-A n y) is zero if x-A fi y is 
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less than zero. 



Since A is a Beta random variable and, 
n ' 



hence, bounded above by one, x-A^y can not be negative if 

x > y. However, if y > x, then A must be restricted to lie 

n 



in the range (0,^) 


. Taking this restriction into account 



and writing the RHS of II.B.3.3 as a convolution 



PIX ni x l X „-l-J" ■ 


L 

= / f A (z)F D (x-yz) dz , (II.B.3.4) 


and 

1 ; 

v y 


if x _> y, 
if x < y. 


where f^(z) is the 


density of a Beta random variable as in 



II. B. 1.3, and F D is the distribution function of a Gamma 
(q,^) random variable. 



Of course, 


to get the conditional density of X 

n 


given X p _ 1 = y, we 


must take the derivative of II.B.3.4 


with respect to x. 


Recognizing that the upper limit may 


be a function of x 


and applying Leibneitz’s Rule where 



appropriate yields 



f | (x) = 

n • n- 1 


L 

/ f, (z) f_(x-yz)dy, (II.B.3.5) 

o A u 


and 

(> 


if x >_ y. 




if x < y. 
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where (x) is the conditional density of X given 

n 1 n-1 n 

Vl ' ^A^ 2 ^ density, and f D (x-yz) is the Gamma 

density as in II. B. 1.2. Writing this result in terms 
of the densities involved we have 




(x) 

n-1 



f Lii— L— z k *3 q n— z ) q ^J^_(y_ vz )3 k(x yz) 

j n r (x-q) F (q) z (i z) r(q) (x yz) e 



dy , 



with the condition on L, as in II.B.3.5. And finally. 



f„ , x (x) = { liM ,)k q e- kx 

A n n-1 T(k-q) [T(q) ] z 



x / z k q 1 (l-z) q 1 (x-yz) q 1 e kyz dy, (II.B.3.6) 
0 

and 

/I if x > y, 

L = 

- if x < y. 

Y 1 



As a check on the derivation of the conditional density, 
the conditional density and conditional expectation were calcu- 
lated for values of k and q which produced simple integrands. 

One of these cases was k = 2, q = 1. Then k-q-1 = 0 and 
q-1 =0. In these parameter values II.B.3.6 reduces to 



f X lx = 2e 2X / e 2yz dz L 

n 1 n- 1 0 



1 



x 

y 



if x >_ y , 
if x < y. 
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After integration. 



J X n-l 



(x) 



-2x 

e 






if x >_ y. 



if x < y. 



(II.B.3.7) 



Since the two expressions in II.B.3.7 are non-negative, we 
can insure that this is a density by verifying that its inte- 
gral is one . 



/ f , v (x) dx 
0 X n |X n-l 



-2x 00 ~ , 2y 

/ [— - — ] dx + / e [ - — ] dx 

o y y y Y Y 



^.1 e 



1 y 1 y _ 2 x e 2y 1 /' -2x 

h / dx “ h / e dx+ " y ] / e ^ 

y o y 0 y y y 



2y 2y 2y y 



= 1, 



as required. We can also take the conditional expectation to 
see if it equals y + j as required by II. B. 3.1. Thus, 



/ xf | (x)dx 

0 X n ' n-l 




-2x 
e 

y 



] dx + / 

y 



xe 



-2x ^e 



2y 

y” 
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GO 



t Xf x | X (x)dx " 
0 n' n-1 



7 /*»* - hf 



y _ 

xe 



2x 



dx + [• 



2y 




xe ^ X dx 

y 




-2y -2y 

e ■* e •* 

2 4 ^ 




4y 




l 

2 ' 



as required. 

Using k = 3 and q = 1 produces a density of 



f X | X < x) = 

n 1 n-1 



2e- 3x [i(e 3 y-S 
y 

2x 2 



3y ± 

w + 3y> 1 if x i y - 



(II.B.3.8) 



2 ^ 2 

y 3y 



/ 1 -3x, 

(1 - e ) 



if x < y . 



This density is non-negative and integrates to one. It also 
produces a conditional expectation of ^ i as required by 
II. B. 3.1. These results are also of use in validating the 
results obtained from numerical integrations of II.B.3.6 in 
estimation applications. 

This conditional density can be used to form the joint 
density of X^ , X 2 , . . . > X^ and, hence, the likelihood function 
of the data. This subject will be addressed in the following 
section . 

4 . Maximum Likelihood Estimation 

Once the conditional density of X n given X n _-, and the 
marginal density of X^ are known, it is possible to evaluate 
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the joint density of X, X . Since the first-order 

-L z n 

autoregressive process is Markovian, as can be seen directly 
rrom the defining equation II. B. 1.1, the following equation 
is valid: 



f (x j I x j-l' x j-2' • * * ' x i ) 



f(x j | x j _ 1 ) 



n > j > 2 



(II .B.4 .1) 



Applying II. B. 4.1 n-1 times to the joint density of X ,X 

n n- 

. . . , produces 



f(x n' x n-l' 



,x) 



f l (x nl 



x n-l )f l (x n-l 



x n-2> 



f 1 (X 2 



x) f 2 (x^) 



(II. B.4 .2) 



where f is the joint density of X n ,...,X^; f^ is the condi- 
tional density of X^ given Xj_^; and f^ is the marginal den- 
sity of the {X^} sequence. 

Viewing the joint density as the likelihood function, 
letting L be the likelihood function, and taking natural 
logarithms of each side of II. B.4. 2 produces 



n-1 

In L = In f 2 ( x 1 ) + _ l In f ] _(x i+1 |x i ) 



(II. 3. 4 .3) 



Recall that in Section II. B. 3 the conditional density 
was determined to be 
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^(x|y) = { 



r(k-q) [T(q) ] 





0 



1 if x >_ y. 



(II. B. 4 .4) 



L 



x . _ 

- if x < y. 



For a given set of data the likelihood can be viewed as a 
function of the parameters k and q. Let 



Assume for the moment that a procedure has been established 
to evaluate G(k,q) given k,q (0 < q _< k) and the data. Consider 
the problem of constructing a program to determine the values 
of k and q which maximize G(k,q) . An outline of the program 
can be constructed as follows: 

1. Read the initial values for k and q. 

Read the data. 

2. Determine a direction of search. Use the following 
difference equations to approximate derivatives. 




(II. B. 4 .5) 




G(k+Ak,q) - G(k,q) 
Ak 



(II.B.4.6) 



fq G < k '<3> 



~ G ( k , q+Aq) - C-(k,q) 
Aq 



(II .B .4 . 7) 
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Let Dk. = Gj. k+ Ak,q) - G (k,q) and = G(k,q+Aq) - G(k,q) 

i Ak M i Aq 

, nv 

for the ith iteration. If we define VG. = (J: i) , then VG. 

1 Dq^ 1 

approximates the gradient of G at the current k,q values. 
For i = 1, let the direction of search, d^, be VG^. For 
i > 1, define the direction of search, d. , to be 



d k - VG i + B k-1 d ]< _ 1 , (II.B.4.3) 

vgTvg. 

where B, is — ~ — , the ratio of the length of the 

l VG i-l 

present and preceding gradients. Formula II.B.4.8 is the key 
equation in implementing the Fletcher-Reeves Conjugate Gradi- 
ent Algorithm. Once d k has been selected, normalize its 
length to one. 

_3 

3. Let the initial step length, SL, be 10 and let N = 1. 
Compute the trial values of k and q, Tk and Tq, using 



Tk = k + SL * Dk^ , 
Tq = a + SL * Dq ± . 



(II .B.4 .9) 



If G (Tk, Tq) > G ( k , q) , let SL = 2 * SL, otherwise go to 4. 
If N = 10; k = Tk, q = Tq, go to 2. (No step larger than 
2 10 * 10~ 3 ~ 1. 0 is allowed.) If n <10; N = N+l, go to 3. 

4. If N > 1 ( at least one step produced an increase), 
use a golden section search between Tk,Tq for step N— 2 and 
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Tk and Tq for step N to determine the maximum function 
value and the k and q values, kMAX and qMAX, which produced 
it. Here k = kMAX, q = qMAX, go to 2. 

If N = 1, go to 5. 

5. Since the initial step along the indicated direction 

produced a decrease in the function value, check to see if 

you are at a local maximum. Determine the function value 

at points at 30° intervals on the circumference of circles 

-3 -2 -1 

with radii of 10 ,10 ,10 (0° is parallel to the q axis) . 

If the maximum of these test values is greater than the pre- 
sent value, set k and q to the values which produced the 
maximum value, set the present function value to the maxi- 
mum, set i = 1 and go to 2. Otherwise terminate. 

The program above assumed that, given q, k and the 
data, the value of the likelihood could be calculated. One 
difficulty in performing the calculation is that the integrand 
of the conditional density may contain singularities. As a 
precondition to using an IMSL routine to evaluate the inte- 
gral, these potential singularities must be removed. The 
technique employed requires that the coefficient of the term 
that goes to infinity as one of the limits of integration is 
approached is added and subtracted. The part that is added 
is then integrated separately and added to the part that is 
evaluated by the IMSL routine. 

To procede with this technique we first split the 
integral into two parts. Thus, ignoring the part of II.B.4.4 
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outside the integral, we have 



/ L z k - < ’- 1 (l- Z ) q ' 1 (x-y 2 ) < 5- 1 e k y z dz = J^z 1 ^- 1 ( l-z)^ 1 

o 0 



x (x-yz) q '*'e^^ Z dz + / z k q ^(l-z) q ^(x-yz) q ^e^^ z dz 

L/2 

(II. 3. 4. 10) 

In the first part of the RHS of II. B. 4. 10 the term Z k -< 5“1 could 

tend to infinity as z tends to zero if k-q-1 <0. If we set 

z equal to zero in the remainder of the integrand we get 

(l_z) < 3- 1 ( x -yz) < 3- 1 e k y z | = x q- ^. Adding and subtacting this 

z=0 

term times the term that contains the singularity, we have 



r L/2 k-q-1,, x q-l, v q-1 kyz, 

J z M (1-z) ^ (x-yz) n e 2 dz 

0 



k-q-1,, ,0-1, , q-1 kyz q-1 k-q-1, q-1 k-q-1, 

= J z M (l-z) M (x-yz)^ 1 e 1 -x n z ^ +x^ z n dz 

0 



L/2 

= / 2 



/2 



k - < 3- 1 [(l-z) < 5- 1 (x-yz) q - 1 e k ^ z -x < 3- 1 ]dz + x^' 1 / z 



k-a-1 



dz 



Thus , 

J L/2 Z k-q-1 ( l-z) q - 1 (x - y2 ) 1-l e k y*dz 
' 0 



= / L/2 z k -'3- 1 ((l-z) q - 1 (x-yz) q - :L e k l' 2 -x q - 1 )dz + x < 3- 1 ^ 

0 



(fo k - q 
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(II. B. 4. 11) 



Recall that since q < k, k-q > 0 and the second integral on 
the RHS can, in fact, be integrated as shown. Now the inte- 
gral on the RHS can be evaluated by IMSL routine DCADRE and 
the second part of the RHS can be easily computed. 

Applying this technique to the second half of the 
integral, recalling that the case where L = 1 and L = — 
must be considered separately produces 

f 1 z k " < J“ 1 (l-z) < J“ 1 (x-yz) q ’ 1 e kyz dz 

1/2 



= f [z k (x-yz) (x-y) q '*'e^ y ] (l-z)^ ^dz (II. B. 4. 12) 

1/2 

+ (x . y) q-i e ^ 



and 

f X/Y 2 k-q-l (1 . 2) q-l (x . yz) q-l e kyz dz 
x/2y 



= J X/Y (x-y Z )^ 1 [ Z k - q - 1 (l-z) q - 1 e kyZ -(2) k - < 3- 1 (l-f) q - 1 e kX ]<3 Z 
x/2y 



(|) q 



(II. B. 4. 13) 



Since q > 0, the second terms in the two previous expressions 
are properly integrated. 
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Two points should be noted. First, if k-q-1 > 0 or 
if <3 - l > 0, then these steps are not required. But whether 
they are required or not, they are always accurate. Since 
the exact path of the search algorithm is unknown at the 
start, these expressions are used throughout to insure accu- 
rate calculations regardless of the k and q values encoun- 
tered. Second, if x = y, two parts of the integrand simul- 
taneously tend to infinity and this procedure breaks down. 

This does not pose a problem for continuous data since the 
probability of this occurring is zero. However, if discrete 
values are used or if the data is truncated because of limits 
on the accuracy of the measuring device, then the data may 
have to be preprocessed to insure that successive values are 
not equal by adding a small increment to one of the values. 

When the program was written, its accuracy was veri- 
fied by three checkes. The case k = q implies independence 
in the basic model since as q tends to k the probability that 
A n equals zero tends to one and the probability that B n 
equals one tends to one. Hence, in the limit X n = G n and G n 
is an iid sequence. The logarithm of the likelihood function 
for independence was calculated and compared to the program 
results for several values of k and q where k = q. The two 
calculations were equal within machine roundoff and compu- 
tational accuracy. The special cases of k = 2, q = 1 and 
k = 3, q = 1 discussed in II. B. 3 were also computed. The 
logarithm of the likelihood function was computed using each 
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of the conditional densities derived in II.B.3.7 and II.B.3.8. 
When the results of these calculations were compared to the 
program results with the specified k and q values, the re- 
sults were equal within calculation and roundoff accuracy. 

The results of the tests of this program when used 
with simulated data with known parameter values are presented 
in Section II ,E . 

Note that there are natural moment estimators for the 
three parameters, k and q and y in this model. These follow 
from the fact that 



P (1) = 

C 2 (X ) 



CORR (X X , ! ) 
n n+1 

Var (X n ) 

[E (X n ) ] 2 




Var (X) 
[E (X) ] 2 




1 

k 



Thus we use for moment estimations 



M = x 



(II .B. 4 . 14) 



q = (l-p(l)k 



(II. B. 4 .15) 



k 




(II. B. 4. 16) 



These moment estimations will be compared to maximum likeli- 
hood estimations in the case where y = 1 in Section II. E. 
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5. Directional Moments 



Unlike processes with normal marginals, non-normal 
processes are not completely determined by their correlation 
structure. Directional moments not only demonstrate this 
difference, but also help to differentiate processes with 
similar correlation structure and identical marginal distri- 
butions. They can also be used to help determine parameter 
values. In addition, they may also be viewed as another 
way of characterizing the joint distribution of the process. 
With 



X 



n 



= A X . + B n G , 
n n-1 n n 



with all random variables defined as in II. B. 1.1, we first 

address E(X X 2 ). We have 
n n-1 



X 



n 



= A X n + 
n n-1 



B G 
n n 



x x 2 . 

n n-1 



A X 3 . 
n n-1 



+ B G X 



n n 



2 

n-1* 



Taking expectations, recalling G^ and Xj are independent if 
i > 1, <'B } and {G } are independent, and A^ and X. are 

independent if i > j . Then 

E(X n X n-l } = E(A n )E(X n-l } + E (B n } E (G n } E (X n-l } 

k-c k(k+l) (k+2) , q t . k (k+1) 

ic 73 T k ,2 
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E(X nVl> 



- >iiIS±lI ( 0sz2l . (k+2) + 



q] 



E (X X , ) 
n n-1 



k(k+l) r k +2(k-q) , 

3 *■ k 

k J K 



(II. B. 5.1) 



Since 



X n = A l x l i + 2 A B G X . + B 2 G 2 , 
n nn-1 nnnn-1 n n' 



2 

X„X„ i 
n n-1 



2 3 2 2 2 

AX, + 2A B G X . + B G X , 
n n-1 n n n n-1 n n n-1 



Taking expectations as before 



E(X n X n-l> ' E<A n ,E(X n-l> + 2E < V E < V E ( V E (X n-!> 

+ E(B 2 )E(G 2 )E(X n _ 1 ) 

(k-q) (k-q+1) k (k+1) (k+2) , 2 (k-q) CT 1 k(k+l) 
k (k+1) k 3 + k ‘k’ 1 ^2 

+ q(q+l) k (k+ 1 ) 1 

k(k-l) k 2 

After simplification we get 

E(X 2 X .) = jk- q) (k+1) (k+2) qk (k+1 ) (II.B.5.2) 

n n-1 .3.3 

k k 

Note that these two directional moments are different, 
indicating that the process is not time reversible. 
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6 - P(x n-H >x n ) 

Another characterization of the joint distribution 
and the time directionality of the process is given by 
p ( x n+ ]_ > X ) . There is a simple analytical solution for 
P (X n+l > X n } in the GLAR (D process. Consider, 



P(X . > X ) = P(A .,X„ + B . G > X ) 

n+1 n n+1 n n+1 n+1 n 



' p(B n + l G „ + l > ll - A n + l )X „> 



(II. B. 6.1) 



Recall that B n+ ^ is Beta(q,k-q) and A n+ ^ is Beta(k-q,q) . 
Hence, [1-A is Beta(q,k-q) . Since G n+1 and X n are inde- 

pendent and have the same marginal distribution and A n+ ^ and 
B n+ ^ are independent, each side of the inequality in II. B. 6.1 
has the same distribution and the random variables are inde- 
pendent. Hence 



P (X 



n+1 



X ) 
n 



0.50. 



This is a strong property of the process. While the 
process, as seen by the directional moments, is not time 
reversible, the fact that p ( x n+ ]_ >x n ) = °- 50 means that 

sample paths will have as many "runs up" as "runs down". 
Sample paths are given in Figures II. B. 6.1 and II.B.6.2. 

An additional point of interest occurs when k = 1 
and the process has exponential marginal distributions. 
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Another exponential first-order autoregressive process in 

which P > X n ) = 0.50 is the PREAR(l) process. This is 

a special case of the two parameter NEAR ( 1 ) exponential 

process of Lawrance and Lewis [Ref. 8 ] in which the two 

parameters a and 6 are related by Q = . The two 

2- a 

exponential processes are very similar in sample path proper- 
ties. However, the GLAR(l) process has a smoother joint 
distribution. In fact, the likelihood for the PREAR(l) 
process is discontinuous, making it difficult to get maximum 
likelihood estimators. 

C. FIRST-ORDER MOVING AVERAGE BETA-GAMMA MODEL, GLMA(l) 

1 . Introduction 

Another special case of the GLARMACp , q) model is the 
first-order moving average model where p = 0 and q = 1. This 
arises naturally from the key result that an {X n > sequence 
can be formed by the random sum of two independent Gamma ran- 
dom variables. In the first-order autoregressive case of 
Section II. B, the generation scheme was given by equation 
II. B. 1.1 and is repeated here 



X 



n 



AX i 
n n-1 



+ B G . 
n n 



The distribution of the (X n ) sequence depends on the inde- 
pendence and distribution characteristics of X n-1 and G n - 
It should be noted, however, that any two independent random 
variables with the required Gamma distribution could be 
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substituted for X n _^ and G n without changing the distribution 

of X . In particular if we substitute G for X , and G , 

11 n n-i n-1 

for G n , we produce the first-order moving average process 
which generates the {X n > sequence using 



X_ 



n 



= A G + B G -i / 
n n n n-1 



(II .C.l . 1) 



where {X n , n = 1,2,...} is a second-order stationary sequence 

of random variables with marginal Gamma(k,l) distributions, 

(A n , n = 1,2,...} is an iid sequence of Beta(k-q,q) random 

variables; {B^, n = 1,2, } is an iid sequence of Beta (q,k-q) 

random variables; (G n , n = 0,1,...} is an iid sequence of 

Gamma (k,l) random variables; (A }, (B }, and (G } are inde- 

n n n 

pendent; 0 < q £ k. The Gamma random variables are para- 
meterized as in II.B.l with density as in II. B. 1.2. The 
Beta random variables have density as in II.B.l. 3. 

In this section we will address the correlation struc- 
ture of the (X } sequences and that of the (X } , (G } se- 
n n n 

quences, theoretical ranges for possible correlations for 

the (X n } sequences, directional moments, and the P (X n+ ^ > X n ) . 

2 . Correlation Structure 

Using II. C. 1.1 to define X and X , we have 
^ n n-i 



X 



X , = (A G + B G -, ) (A^ i G , + B ,G _) 

n n-1 n n n n-l n-1 n-I n-1 n-z 



= AA n G G„ t + A B ,GG - + A , B G ,+BB , G ,G„ , 
n n-1 n n-1 n n-1 n n-2 n-1 n n-1 n n-1 n-1 n-2 
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Using the iid nature and independence of {A } , {B }, and 

n n 

{ G n ) an< ^ taking expectations 

E (X n X n-l ) = E (A n } E (A n-l } E ( V E (G n-l } +E (A n } E ^n-l 5 E (G n } E (G n-2 } 

+ E (A n-l> E (B n> E (G n-l» +E (B n> E < B n-l> E < G n-l> E <G n-2> 

= (^) 2 + (!^)(|) + (SjS) (2) ( ktk+il) + ( a, 2 

k 

= 1+ ( k } (II. C. 2.1) 

Therefore, 

C0V(X n' X n-l> " 

and 



CORR(X n ,X n _ 1 ) = (1-|)(2). (II.C.2.2) 

A simple calculation will show that for lags greater 
than one the correlation is zero. So equation II.C.2.2 plus 
the knowledge that greater lags are zero is sufficient to 
specify the correlation structure of the {X^} sequence. 

One might note at this juncture that the correlation 
of the (X n > sequence is constrained to lie in the interval 
(0,j). The reason for this constraint and a method for re- 
laxing it will be discussed later. It is also worthy of 
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note that the {X n > sequence is stationary and has the same 

marginal distribution as the {G n > sequence, if X Q = G Q . 

As indicated in II. B. 2 the {X } and {G } sequences 

n n 

can be considered to be a bivariate, correlated Gamma(k,l) 
process. As such, the correlation structure of this bivari- 
ate Gamma may be of interest. Consequently, we first develop 

the correlation of X and G in the standard fashion. 

n n 



X 



n 



AG + B G -I 
n n n n-1 



X G _ 
n n 



AG + B G G . 
n n n n n-1 



Taking expectations as before. 



E(X n G n ) - E(A n )E<G„) + E <V E(G n )E(G n-l> 

= (^3.) ( klk 2 ~-) + | 

k 2 + k - q 
k 2 

Therefore, 



C ° V <V G n> 




and 
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(II .C.2 . 3) 



CORR(X n ,G n ) = ^3. 



Now consider the correlation of X and G , . We have 

n n-l 



X„ = AG + B G . 
n n n n n-l 



X G , 
n n-l 



A G G , + B G , 
n n n-l n n-l 



Expectations in the standard fashion produce 



E(x n G n-l> 



E (A ) E (G ) E (G ,) + E(B >E(G* ,) 
n n n-l n n-l 



_ k-q q , k [k+1] x 

" k k ( k 2 > 



k +q 
2 



2 

k‘ 



Hence , 



C °V (x n ' ^n-1 5 = 



_S_ 

k 2 



And finally 



C0ERt V G n-l> = k 



(II. C.2. 4) 



A simple calculation convinces one that the correlation for 
lags greater than one is zero. In addition, it is clear 
that CORR(X n ,G n+m ) = 0 for m = 1,2,... Hence, II. C.2. 3 and 
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II.C.2.4 are sufficient to specify the correlation structure 
of the {X n } and { G n } sequences. 

It has been noted before that the range of correla- 
tions for the first-order moving average process generated 
by the Beta-Gamma method of II. C. 1.1 is constrained to lie 
in the interval from zero to one-quarter. There may be other 
random linear combinations of Gamma random variables which 
give a moving average process with Gamma marginals and a 
greater range of positive correlations. Thus we now examine 
a more general hypothetical generation process to prove that 
any random, linear combination of two independent Gamma random 
variables which generates a sequence with the same first two 
moments as those Gamma variables has a correlation that lies 
in this same interval. In fact, this proof only requires 
that the dependent random variable have the same first two 
marginal moments as the generating Gamma variables. 

THEOREM: 

If the {X^} sequence is generated by 

X = A G + B G , , (II.C.2.5) 

n n n n n-l 

where {X } is a second-order stationary, non-negative sequence 
n 

of random variables with the same first two moments as the 

{G } sequence; {A } and {B } are iid sequences of random 
n n n 

variables; {G } is an iid sequence of Gamma (k,l) random 
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variables; and (A n ), ( B n K and (G n ) are independent, then 

0 < CORK ( X , X ,) <0.25. 

— n n-l — 

Proof: Since {A_}, {B_}, and (G } are independent and {X } 

n n n c n 

is non-negative, {A } and {B } must be non-negative. Hence 

n n 

E (A) >0 and E (B) > 0. 



Hence, 0 <_ E(A) <_ 1 and 0 <_ E(B) £ 1. 

Computing the serial correlation of fX n ) yields 



Taking expectations of II.C.2.5, we have 



E(X ) 
n 



E(A)E(G) + E (B ) E ( G , 
n n n n— l 



1 = E (A) + E(B) 



(II.C.2.6) 





Then 



E (X X , 
n n-l 





1 + E (A) E (B) (i). 
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1 , 



Since E(X ) = E (G ) = 

n n 

COV(X n ,X n _ 1 ) = E(A)E(B)(i) 

And since VAR(X ) = VAR(G ) , 

n n 

CORR(X n ,X n _ 1 ) = E (A) E (B) 

Using II.C.2.6 and its consequences, we have 

0 £ CORR(X n ,X n _ 1 ) = E (A) [1 - E (A) ] < j. (II.C.2.7) 

So, in general, if {X n } is second-order stationary with the 
same first two moments as {G^}, the serial correlation of 
{X^} is bounded below by zero and above by one-fourth. Q.E.D. 

This constraint on the correlation appears to be 
restrictive since, in the classical case, when two normally 
distributed random variables are added to produce a normal se- 
quence, the range of correlations is (- ^ ,^) . The two situa- 
tions, however, are not comparable. It is clear upon reflection 
that the constraints imposed on the (X n ) sequence in the pre- 
vious theorem are more severe than those imposed upon the 
classical normal case. In the above theorem we required that 
both the mean and variance of the { X n } sequence equal that of 
the innovative sequence. However, in the classical normal 
case (where zero mean normals are used as innovative factors) 
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only the mean of the generated sequence is equal to the mean 
of the innovative sequence. The variances are not equal. 

We now examine the case where the generated sequence is re- 
quired to have the same mean as the innovative sequence, but 
is free to have a different variance. (This is the case 
with the usual constant coefficient, linear additive MA(1) 
scheme . ) 



II.C.2.5 and all variables are defined as for that equation 
except that (X n ) is only constrained to have the same first 
moment as -[ G n } , then 0 £ CORE (X n , X n _^) £ 0.5. 

Proof : Taking expectations of II.C.2.5 with the new circum- 
stances produces 



and 0 < E (A) £ If 0 £ E (B) £ 1 by following reasoning identi- 
cal to that above. Calculation to determine the serial corre- 
lation can initially proceed as usual. 



THEOREM 



If the non-negative {X n > sequence is generated by 



E(X) 



[E (A) + E(B)]E(G) 



1 = E (A) + E (B) 






1 + E (A) E (B) (jjr) 
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C0V(X n' X n-l> 



= E (A) E (B) (£)• 



To this point all calculations and reasoning have been the 

same as that which produced II.C.2.7. However, since {X n } 

is not constrained to have the same second moment as {G } 

n 

the most explicit result obtainable is 

E ( A) E (B) (i) 

C 0 RE( V X n-l> “ V A R DO ' (II.C.2.8) 



where VAR(X) is, of course, a function of VAR(A) , VAR(B) , 

and VAR(G) . Since it is obvious that the smaller the value 

for VAR(X) , the greater the serial correlation for {X }, 

n 

let us reduce VAR(X) to its smallest values. Since the dis- 
tribution of {G^} has been specified, its variance is fixed. 
Let P (A n = a) = 1 and P(B n = b) = 1. Then trivially E (A) = a, 
E (B) = b and VAR(X) = (a 2 + b 2 ) (j^) . Under these conditions 
II.C.2.8 becomes 



CORR(X n ,X n _ 1 ) 



2 

a 



ab 



+ b 



2 * 



If we further specify that a = b, then CORR(X n ,X^_^) achieves 

its maximum of 1/2. Q.E.D. 

The situation developed above is comparable to the 

classical situation where the innovative factors have distri- 
2 

bution, N(0,a ). Except for the degenerate case where one 
coefficient is zero, the sequence generated by a linear com- 
bination of innovative factors may have the same first moment 
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as the innovative factors, but will have a different second 

moment. So, under comparable conditions the random, linear 

combination of Gamma(k,l) random variables can produce 

positive correlations equal in magnitude to the positive 

correlation produced by a classical normal process. The 

distribution of the {X } under these circumstances is unknown. 

n 

If the distribution of the {X } is constrained to 

n 

be Gamma (k,l) and the Beta-Gamma generation scheme is used 
to generate the {X n }, then the maximum correlation that can 
be achieved is one- fourth. 

3 . Directional Moments 

As mentioned in II. B. 5 the directional moments of a 

non-normal process are not necessarily equal and may provide 

valuable information about a time series. First, consider 

E (X 2 X . ) . From II.C.l. 1 
n n-l 



2 2 2 2 2 
X„ = A„G„ + 2A B G G , + B^G^ . 
n n n n n n n-l n n-l 



and 



X = A . G i+B . G 0 

n-l n-l n-l n-l n-2 



Therefore, 



X 2 X . 
n n-l 



2„3 



A n A n- l G n G n-l + 2A n A n-l B n G n G n- 1 +A n- l B n G n-l +A n 3 n- 1 ’n G n- 2 



2 2 

+ 2A B B G G i G ~“t~B B G G ^ 
n n n-l n n-l n-2 n n-l n-l n-2 
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Taking expectations as in II. C. 2 yields 



E(X 2 X ) = (k ~ q) (k-q+1) + 2 ( k-q ) _ (k+1 ) + (k-q) q (q+1) (k+2) 

n n_1 k 3 k 3 k 4 



, (k-q) (k-q+1) q . 2(k-q) 2 q 2 , q 2 (q+1) 

v 3 + ,3 ,3 

k k k 



Upon simplification this produces 



E(X X . . ) 
n n+1 



= Ijl^tlk-q+Sl+i ^a t ^k+ll+k+l j 



+ ^r[2k-q+l] 
k 



(II .C. 3.1) 



In an analogous fashion 



XX 2 . = A A 2 . G G 2 . + 2A A .B . G G ,G n , + A B 2 , G G 2 - 
n n-1 n n-1 n n-1 n n-1 n-1 n n-1 n-2 n n-1 n n-2 



+ A 2 . B G 3 , + 2A .B B ,G 2 ,G n , + B B 2 i G , G 2 - 
n-1 n n-1 n-1 n n-1 n-1 n-2 n n-1 n-1 n-2 



Taking expectations we have 



E (X X . ) 
n n-1 



(k-q) 2 (k-q+1) , 2 (k-q) 2 q , (k-q)q(q+lj_ 

1 + 3 3 

k 4 k^ k J 



(k-q) (k-q+2)q(k+2) , 2 (k-q) q 2 (k+1) , q 2 (q+1) 

+ k 4 k 3 ' 



which simplifies to 
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E ( XX ,) 

n n-1 



(k-q) 3 (k+q+1) + q (q+1) ( 2k-q) 
k 3 k 3 

+ (k-q) qjk+1) 2 



(II.C.3.2) 



4. Empirical P(X , . > X ) 

n+1 n 

No analytical procedure was found to determine 
P(X n+1 > X ) . Hence, a simple computer program was constructed 
to evaluate this condition for a series of sixty-eight thou- 
sand pairs of numbers generated by the Beta-Gamma scheme for 
each of ten random number seeds. The answer obtained was 
considered to be accurate within 0.001. The comparisons were 
run for each of seventy-nine values of q, from 0.05 to 3.25 
in increments of 0.05. All of the results of the comparisons 
fell in the range 0.499 to 0.501. Fourteen of the values were 
different from 0.500. No pattern was apparent in the devia- 
tions from 0.500 and these deviations were considered to be 
random fluctuations within the given margin for error. It 
thus seems clear that P(X n+ ^>X n ) for this process is like 
the GLAR(l) process but no proof has been found. 

D. OTHER CASES OF THE GLARMA ( p , q ) MODEL 
1 . Introduction 

A primary advantage of the GLARMA (p,q) model is the 
ease with which it can be adopted to cover a variety of 
special cases. Two special cases, the first-order autoregres- 
sive GLAR(l) and the first-order moving average GLMA(l), were 
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covered in II. B and II. C. The intention here is to briefly 
present three additional cases of the general model and derive 
the correlation structure of each case. The special cases 
considered are the first-order mixed model, GARMA(1,1), 
the second order autoregressive model GAR (2) , and a bivariate, 
first-order, autoregressive model BGAR(l) . The purpose 
in presenting these cases is to demonstrate the flexibility 
of GLARMA(p,q) and not to present a complete, detailed dis- 
cussion of each model. Further extensions of the special 
cases of the GLARMA(p,q) model from the examples given are 
obvious. Details will not be given. 

2. GLARMA (1,1) 

Consider the following scheme for generating an {X^} 
sequence of random variables. 



X 



n 



— BA -,+CG, 
n n-1 n n 



n 



= DA , + F G . 
n n-1 n n 



(II. D. 2.1) 
(II .D.2.2) 



where {X , n = 1,2,...} is a second-order stationary sequence 
n 

of Gamma (k,l) random variables; {A^, n = 0,1,...} is a 
second-order stationary sequence of Gamma (k,l) random 
variables; (B^, n = 1,2,...} is an iid sequence of Beta(k-q,q) 
random variables; {C^, n = 1,2,...} is an iid sequence of 
Beta(gfk - cj) random variables; (D^, n = 1,2,...} is an iid se- 
quence of Beta(k-r,r) random variables; {F n , n = 1,2,...} is 
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an iid sequence of 



Beta(r,k-r) random variables; {B }, (C }, 

n n 

(D n ), {F n ), and (G n ) are mutually independent; 0 < q £ k; 

0 < r £ k. The Beta(m,n) density is given in II. B. 1.3; the 

Gamma(k,l) density in II. B. 1.2. 

Before the serial correlation of the {X^} sequence 

can be determined, the serial correlation structure of the {A } 

n 

sequence and the correlation structure between the {A^} and 
{G n } sequences must be derived. Proceding first with the 
serial correlation of the {A n > sequence, from II.D.2.2 we 
have 



n 



= DA . + F G 
n n-1 n n 



So 



A A . 
n n-1 



DA . + F G A , 
n n-1 n n n-1 



Using the iid nature of {D n > , {F n >, { } ; noticing that when 
i > j, D i and Aj , F^ and Aj , and G^ and Aj are independent; 
recalling the independence of (F n ), {G n }; anc ^ taking expec- 
tations yields 



E < A nVl> 



,k-r N k (k+1 ) , r 

( k ‘ k 2 k ' 



E (A A 
n 



n-l J 



k +k-r 



(II .D.2 .3) 
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Therefore , 



C 0 V(A n-Vl> 



k-r 



and 



CORR(A n ,A n _ 1 ) = (II.D.2.4) 

In fact A n is just a GLAR(l) process, so the result is not 
surprising. Using II.D.2.2, II.D.2.3, and an induction argu- 
ment leads to the general m-step correlation formula 

CORR ( A , A ) = n > m > 0. (II.D.2.5) 

n n-m x. 

The correlation structure between {A } and (G } can 

n n 

be derived in a similar fashion. However, it is more direct 
to note that since the {A n > sequence is the same as the GLAR(l) 
process of Section II .B, the {A n > , {G n > correlation structure 
will be the same as that derived for the {X n >, {G^} sequences 

in II. B. 2. Hence, 

CORR(A n ,G n _ m ) = (£) n>m>0. (II.D.2.6) 

Of course, if j > i, then CORR(A^,Gj) = 0 

Now the serial correlation for the {X n > sequence can 
be found. From II. D. 2.1 
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X n = 


BA , + C G 
n n-1 n n 


X n X n-l = 


(B A i + C G ) (B , A 9 + C„ ,G„ n ) 
n n-1 n n n-1 n-2 n-1 n-1 


= 


B B , A t A 0 +B C A ,G n +B .C 0 G 

n n-1 n-1 n-2 n n-1 n-1 n-1 n-1 n n-2 n 

+ C C , G G , . 
n n-1 n n-1 



Using the stationarity of the {A n > sequence, the iid nature 

and independence of {B }, {C } , {G }, the fact that A. is 

n n n 1 

independent of Gj when j > i, and the fact that {A n } is inde- 
pendent of (B n ) by construction and taking expectations, we 



have 




E ( XX ,) 
n n— l 


■ EIB ») EIB n-l )E( VlV2 )tEIB n IEIC n-l 1E( VlVl l 

+ E < B n -l )E<C n )E(A n-2 )E<G n )+E(C n )E(C n-l >E(G n )E(G n-l ) ' 
= [E(B n )] 2 E(A n . 1 A n . 2 ) + E(B n )E(C n . 1 )E(A n . l G n . 1 ) 

+ E(B , ) E (C ) E (A -,)E(G) + [E(C) ] 2 [E(G) ] 2 . 
n— 1 n n-l n n n 

(II .D.2.7) 



From II.D.2.1 

E(X n )E(X n _ 1 ) = [E(B n )E(A n _ 1 )+E(C n )E(G n ) ] x 

[E (B . ) E (A ~ ) +E ( C , ) E (G ,)], 
n— l n— z n— ± n — - l 
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E(X )EU ,) 
n n-l 



= [E(B ) ] 2 E (A , ) E (A p) + 2E (B ) E (C ) E (A ) E (G ) 
n n - ± n— z n n n n 

+ [E(C n ) ] 2 [E(G n ) ] 2 . (II.D.2.8) 

Using II.D.2.7 and II.D.2.8 to compute COV(X ,X .) yields 

n n-1 

C0V «n'Vl' - [E(B n>l 2 < E <\-lV2>- E(A n-l )E < A n-2 )I 

+ [E (B n ) E(C n ) ] [E (^ n _ 1 ^ n _ 1 ) -E (A n ) E (G n ) ] 

Since (A } , {G } , and { X_ } are all Gamma(k,l), VAR(X ) = VAR (A ) 
n n n n n 

= VAR(G ) . Hence, 
n 

CORR(X n ,X n _ 1 ) = [E(B n ) ] 2 CORR(A n ,A n _ 1 )+[E(B n )E(C n ) ]CORR(A n ,G n ) 

From II.D.2.4 and II.D.2.6 we know this equals 

CORR(X n ,X n _ 1 ) = (^) 2 (^) + (^) (|) (f) (II.D.2.9) 

Using II.D.2.5 and II.D.2.6 in an induction argument 
yields the general m-step correlation of 

C0RB(X »' X n-m) - (^)”- 1 1 (^» 2 (^» + (^,(a,<£, ] . 

Recognizing the expression in brackets as CORR(X n , X and 

letting p equal this correlation we have 
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CORK (X ,X ) 
n n-m 



(II. D. 2. 10) 



,k-r,m-l ^ ^ . 

= (-£— ) p , n _> m _> 1. 

Figure II. D. 2.1 shows the possible combinations of 
one- and two-step correlations for the GLARMA (1,1) model. 

This concludes the development of the correlation structure 
of the GLARMA (1,1) model. 

3 . GLAR ( 2 ) 

Jacobs and Lewis [Ref. 12] first developed the follow- 
ing mixture scheme for generating a p-order autoregressive 
processes. We now adopt that scheme for generating a second- 
order autoregressive sequence of random variables. This is 
the special case of GLARMA (p,q) with p = 2 and q = 0. As such 
it closely resembles the GLAR(l) process. Let 



X = B X _ + C G , 
n n n-T n n n 



(II .D. 3.1) 



where {X n , n = -1,0,1,...} is assumed to be a second-order 
stationary sequence of Gamma (k,l) random variables; {B n , 
n = 1,2,...} is an iid sequence of Beta(k-q,q) random varia- 
bles; {C n } is an iid sequence of Beta(q,k-q) random variables; 
{G n } is an iid sequence of Gamma (k,l) random variables; {B^} , 
{C n }, {G n } are independent; also T n is iid with ?{T n =l) = 
l-P(T n =2) = a. The Gamma (k,l) and Beta(m,n) densities are 
found in II. B. 1.2 and II.B.1.3, respectively. 

This generation scheme works even though X n _^ and 
^ are dependent random variables. The mixture of X^_^ and 
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X n _2 produced by II. D. 3.1 is Gamma distributed and indepen- 
dent of G n . Hence, II. D. 3.1 is simply another example of the 
random sum of two independent Gamma random variables producing 
another Gamma random variable. 

Two special cases of II. D. 3.1 are as follows. When 
a = 1, the GLAR(2) process reduces to the GLAR(l) . When q = k, 
the (X n > sequence is iid. 

The serial correlation of the {X n } sequence can be 
calculated in the usual fashion, assuming stationarity of 
the process we have 



X n ’ B n X n-T + C n G n 
n 



k>0; 0 < q < k 



and 



XX . = BX X . + C G X .. 

n n-1 n n-T n-1 n n n-1 

n 



Using the independence of {C n } and {G n }, the fact that X^ 
is independent of C ^ , Gj and B j when j > i and taking expec- 
tations we have 



E (XX ,) 
n n— i 



= aE(B n )E(X^_ 1 ) + (l-a)E(B n )E(X n _ 1 X n _ 2 )+E(C n IE(G n )E(X n _ 1 ) 



- C (tltiil) + (!-„) ( x =3 )E (x n _ 1 x n _ 2 ) + |, 



since E(X n ) = E(G n ) = 1 by assumption. Using the second-order 
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stationarity of { X^} 



E ( X X , ) (1 
n n-1 



(1-g) (k-g) , 
k 1 



a (k-q) (k+1) q 

k 2 k 



Upon simplification 



E ( X X , ) 
n n-1 




+ak-gkq-aq+kq 
k (q+gk-gq) 



Hence, 



COV ( X , X ,) 
n n-1 



(i.) r a (k-q) , 
V q+g (k-q) J 



and 



CORK (X ,X ,) 
n ' n-1 



g (k-q) 
q+g (k-q) 



(II.D.3.2) 



If g = 1, this equals 1 - if k = q then it is zero. Since 
q > k it is clearly non-negative. 

The lag two correlation can be calculated in a 
similar fashion. We have 



X = B X m + C G ; 
n n n-T n n n 



XX _ = BX m X _ + C G X . 

n n-2 n n-T n-2 n n n-2 



and 



71 



E(X n X n-2 } " 



a(^2)E(X„ , X 0 ) + (1-a) (^#)E(X 2 ,) 



k-c 



n-1 n-2 



n-2' 






Using the second-order stationarity of the { X^} sequence we 
can write 

E(X n )E(X n-2> “ “ ( ^> E < X n-l )Et V2> + <1 ^ )< ^ ,|ElX n-2 )l2 

+ k [E(X n )l2 ' 

so that 

C0V(X n' X n-2» = “‘T 2 ’ [E(X n-l X n-2>- E(X n-l )E(X n-2 n 

+ (1-a) (S^2) [E(x^_ 2 )-{E(X n _ 2 )) 2 ] 

^-i‘ E <VJ 2 

- o( T 2 )' E ( X n-l X n-2 ) - E(X „-l ,E ( X „-2) 1 

+ ( ^ 1- ^ -- 2-L [ e (X^_ 2 ) - (E(X n _ 2 )} 2 ] 

= a(^)COV(X n . 1 ,X n . 2 ) + <i=^2ivAE(X n . 2 ) . 



Hence , 



CORR(X n ,X n _ 2 ) = a(n^ 2 )CORR(X n _ 1/ X n _ 2 ) + 



(1-a) (k-q) 
k 
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C°RE(X n ,X n _ 2 ) 



( £l2) , q+ak-2aq . 
' k J l q+a (k-q) J 



(II.D.3.3) 



The solution procedure for II.D.3.3, if followed for 
X n-m' w -“-H Produce the general recursion equation (Yule-Walker) 
that can be used along with II.D.3.2 and II.D.3.3 to compute 
the m-step correlation. The formula thus produced is 

C°RK(X n ,X n . m ) = (^) [o.CORE(X n ,X n+1 _ m ) 

+ (l-ct)CORR(X ,X J ] , (II.E.3.4) 

n n+2-m 

n >_ m >_ 1 . 

As mentioned in previous sections the (X ,G ) pair 

n n 

can be considered to be a correlated, bivariate pair of 

Gamma (k,l) random variables, Therefore, we proceed to 
derive the correlation structure between these two sequences. 
From II. D. 3.1 we have 



X = 3 X m + C G 

n n n-T n n 

n 



and 



XG = 3 X _G + CG. 

n n n n-T n n n 

n 



Recalling the independence of {C^} and {G n > and X^ and G ^ 
when j > i, taking expectations yields 
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E (X G ) 
n n 



,*=£) + (3.) (*> Ogl) 



k 2 +q 

k 2 



Hence , 



and 



COV(X ,G ) = 



n n 



CORE(X n ,G n ) = 2 



(II.D.3.5) 



and 



Continuing this process we have 



= B X + C G 
n n n-T n n 



XG , = BX m G ,+CGG , 

n n-1 n n-T n-1 n n n-1 



Thus 



E(X n G 



n-1' 






(1-cO (^) 



+ 



q 

k' 



Hence , 
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And 



C0V( V G n-l> - a GG ^ ^ X n-1 ' G n- 1 ^ 



CORR(X n ,G n _ 1 ) = a(^)(|). 



One further step in this process using an arbitrary 
m-step lag produces the general recursion formula 

C0RE(X n' G n-m> = <^HaCORR(X n ,G n+1 _ m > 

+ (l-a)CORR(X ,G ^ )], n > m > 2 (II.D.3.7) 

n n+2-m — — 

Figure II. D. 3.1 displays a plot of the possible com- 
binations of CORR(X ,X ,) and CORR(X ,X -,) . Note that when 

n' n-1 n n-2 

a = 1, the GLAR(2) process reduces to the GLAR(l) and 

v- a o 

CORR(X ,X _) = ( — - 1 ) which defines the lower boundary of the 
n n— 2 K 

plot. When a = 0, CORR(X ,X .) = 0 and CORR(X ,X„ 9 ) = 

n n-1 n n-z is. 

which goes front zero to one. In the interior of the graph, 
when a does not assume an extreme value, CORR(X n ,X^_ 2 ) does 
not reach a value of one. This is demonstrated by the 
following calculation. From II.D.3.3 



CORR (X , X ) 
n n-z 



/ k-q ^ [ . q+ak-2gq 1 
' > ' q+a (k-q) J 



If this correlation is to equal one, then 
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(II .0.3.8) 



q+ak-2ctq _ k 
q+a (k-q) k-q 

which quickly reduces to 

q (-2ak-q+2kq) = 0 

Since 0 < q, this requires 

-2ak-q+2aq = 0 



Therefore, 



a = 2 

2 (q-k) 

Since we know that 0 < q £ k, a must be less than zero. 

However, a is a probability and, hence, is non-negative. 
Therefore, the original requirement in II.D.3.8 cannot be 
satisfied. Hence, CORR (X R ,X n _ 2 ) cannot equal one. 

4 . BGAR(l), Bivariate Model 

To this point the only examples of bivariate Gamma 
processes presented were those in which the innovation se- 
quence was one part of the bivariate process with the gener- 
ated sequence the other part. The simple random, linear 
structure of the GLAR(l) process makes it easily extendable 
to a variety of bivariate models. We address only the 
simplest. Consider the following pair of random variables, 

both of which are formed from the same innovation process, (G }: 

n 
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(II .D.4 .1) 



X 




n 



Y 




(II .D.4 .2) 



n 



where (X^, n = 0,1,...} is a second-order stationary sequence 
of Gamma (k,l) random variables; {Y n , n = 0,1,...} is a 
second-order stationary sequence of Gamma (k,l) random varia- 
bles; (B^, n = 1,2,...} is an iid sequence of Beta(k-q,q) 
random variables; {C n , n= 1,2,...} is an iid sequence of 

Beta (q, k-q) random variables; (D n , n = 1,2, } is an iid 

sequence of Beta(k-r,r) random variables; (F n , n = 1,2,...} 

is an iid sequence of Beta(r,k-r) random variable; {G , 

n 

n = 1,2,...} is an iid sequence of Gamma (k,l) random 

variables; (B } , (C }, {D }, {F }, and (G } are independent; 
n n n n n 

0<r<_k; 0 < q <_ k. 



a more general case the {X n } and {Y n } sequences could have 

separate, correlated innovation sequences instead of sharing 

a single {G } sequence. In addition, the {B } and (D } 
n n n 

sequences and the {C n } and {F n } sequences could be correlated. 



{Y n }, it will be necessary to address the correlation of each 
of these sequences with the {G n } sequence. This is most 
easily handled by recognizing that the relationship of the 
{X n } and (Y n } sequences to the {G n } sequence is exactly the 
same as the {A n } sequence to the {G^} sequence in II.D.2.2. 



This is a special case of a general situation. In 



Before examining the correlation between {X n } and 
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Hence, the correlation structure will be the same. There- 
fore, if we let p = CORK (X , X ,) and p v = CORR(Y ,Y .), 

X n n-1 Y n n-1 

n = 1,2,..., these correlation structures can be written 
without further analysis as 



CORR(X n ,G n ) 

CORR (X , G -,) 
n n-1 

CORR (X , G ) 

n n-m 



£ = 1 - CORR(X n ,X n _ 1 ) = 1 - p x ; (II.D.4.3) 

= (g) (^) = (1-P X )P X ; (II.D.4.4) 

= (^) = (1 - P x ) P x , n :> m :> 0 . (II.D.4.5) 



and 



CORR(Y n ,G n ) = 
CORR(Y n ,G n _ 1 ) 

C 0 RE( V G n-m> 



r 

k 



1 - CORR ( Y , Y ,) = 1 - p„; 

n' n-1 Y 



,r . ,k-r. , . , 

k k ~ ^ ~ Py^ Py ' 

= = (1- p y )p” n > m > 0 



(II. D. 4 .6) 
(II.D.4.7) 
(II. D. 4 .8) 



The assumption is that the bivariate process is stationary, 
although it should be noted that starting the univariate 
processes in a stationary mode does not make the bivariate 
process stationary. The initial pair {Xg,Yg} must have the 
bivariate Gamma distribution associated with the stationary 
Markovian process. 

Now we address the cross correlation between (X n ) and 
(Y n ). Start with II. D. 4.1. We have 
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X = B X , + C G 
n n n-1 n n 



and 



XY . = BX , Y ,+CGY . . 

n n-1 n n-1 n-1 n n n-1 



Using the independence of (C n > and {G n } and that of Y 
Gj and Cj when j > i and taking expectations 



E (X Y .) = E (B ) E (X , Y .) + E(C )E(G )E(Y , 

n n-1 n n-1 n-1 n n n-1 



so that 



E (X Y , ) 
n n-1 



- + !■ 



Now starting with II.E.4.2, we have 



Y = D Y , + F G 
n n n-1 n n 



and 



XY = DY .X + F G X . 
n n n n-1 n n n n 



Taking expectations as above 



E(X n Y n ) = !IVEIVn-l' + E(F n )E(G n X n ) 



. and 
l 



) , 



(II .D.4 .9) 
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so that 



Using II.D.4.3 we deduce E (G X ) = — ^3- 

n n ,2 

k 



E(X n Y n ) = + 

k 



(II .0.4 



Invoking the second-order stationarity of the {X n ,Y n } se- 
quences, using II.D.4.9 and II. D. 4. 10 and letting w = E(X i Y i ) 
and z = E(X^Y^_^), we have the two equations 



w . <^>z + <f,<*-±2>, 

k 



(II. D. 4 



z . (^»w + 2. 



(II .D . 4 



Using II. D. 4. 12 to substitute for z in II. E. 4. 11 yields 



w = (fe£)[(^)w + |] + (|)(t+a) 

k 

= [ k 2 -kr-kq+gr 3w + (k-r)g + ( kj+g. ) r 
k 2 k 2 k 3 



When E (X Y ) is substituted for w after simplification, this 
n n 

produces 



<Vn> 



2 2 
_ k q-kqr+k r+rq 

~ k(kq+kr-qr) 



(II. 0.4 



Hence 



COV (X , Y ) 
n n 



ra 



k (kq+kr-qr) 



. 10 ) 

. 11 ) 

. 12 ) 



.13) 
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and 



CORR(X n ,Y n ) 



£3 

kq+kr-qr ' 



k>0, 0 < q _< k, 0 < r £ k. 

(II. D. 4 .14) 



d-P y ) (1-P Y > 
(1-P X )+(1 “ P Y ) P X 



(1-P X ) (l-p y ) 
(l-p y )+(l-p x ) p y ' 



(II .D. 4 .15) 



This latter expression follows since for a given k the corre- 
lation structure is parameterized by r and q or equivalently 
by the serial correlations p and p given at II.D.4.3 and 

X. x 

II. D. 4 .6 . 

We can now substitute II. E. 4. 13 into D.E.4.12 to 

solve for E (X Y . ) . 

n n— l 



E(X Y 1> 
n n— 1 



(*12) r k 2 q-kgr+k 2 r+rg 1 2 

v k ; 1 k (kq+kr-qr) J k 



k 2 q+k 



3 2 2 

r+kqr-k qr-q r 

(kq+kr-qr) 



Hence, 



cov(x n . 



Y n-1> 



. 2 
k q r ~q r 

(kq+kr-qr) 



and 



CORR(X , Y ) 
n n-l 



( 2£ 

kq+kr-qr 




CORR(X Y ) P (II. D. 4. 16) 

n ii x 
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Continuing in this vein produces the general formula 



CORR(X n 



Y 

n-m 



( 2£ ) ( j£ Z3) m 

v kq+kr-qr' v k 



(II .D.4 .17) 



= CORR(X n/ Y n ) p™, n > m > 0. 



Solving for correlations where the {X } lags the 

n 

{Y n ) sequence is similar to the above process, but somewhat 
abbreviated. Starting with II. D.4. 2, 



n 



= D Y . 
n n-1 



F G , 
n n 



we have 



Y n x „-1 



D Y .X i + F G X . 
n n-1 n-1 n n n-1 



Taking expectations as before gives 



Et Vn-l> 



S(D n )D(Y n ,X_ ,) 
n n-i n-l 



E(F n )E(G n )E(X n . 1 ) 



Using the second-order stationarity of (X n ) and {Y n >, 

E (Y .X . ) is known from II. E. 4. 13, so 
n-1 n-1 



E ( Y X , ) 
n n-1 



2 2 

, k-r > , k q-kqr+k r+rq . r 

k k[kq+kr-qr] ] k 



3 3 2 2 

k q+k r+kqr-k qr-qr 

k^ (kq+kr-qr) 
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Hence 




. 2 
kqr-qr 



k 2 (kq+kr-qr ) 



and 



CORK ( Y , X . 

n n-1 



kq+kr-qr ; v k 







Further computations of this nature produce the 
general formula 



To examine the correlation further, note that if 
p x = P y = P/ then II. D. 4. 15 yields 



Thus if p = 0 (i.e., the {X^} and {Y n } processes are iid 
sequences) , this correlation is one. For { } and (Y n ) to 
be iid, we must have X = G = Y . Therefore, this limiting 
correlation does make sense and suggests that perhaps a bi- 
variate Gamma should be used as the innovation process. This 



CORR ( Y , x 

n n— m 



2£ ) ( *z£) m 

kq+kr-qr ; v k ; 



CORRtX^Y^ p™ n>_m>_0; k> 0; 0 <q<_k; 

0 < r < k. 



CORR(X n ,Y n ) 



1 - P 

1 + P* 
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would allow for separate control of the serial correlation 

(auto-correlation) and cross-correlation of the {X } and {Y } 

n n 

sequences . 

If p -*■ 1 (i.e., X n ~ X n-1 , Y n ~ Y n-1 ) , the effect of 

the innovation sequence is slight and the cross-correlation 

between {X n > and {Y^} goes to zero. In a more complicated 

model than we have addressed here the cross-correlation may 

be controlled by imposing a correlation on the {B } and {D } 

n n 

sequences . 

Cross-coupled processes, as discussed in Gaver and 
Lewis [Ref. 2], are possible. These can be used to create 
negative serial correlations in the {X n > and { } processes. 

E. NUMERICAL CONVERGENCE OF THE MAXIMUM LIKELIHOOD COMPUTER 
PROGRAM AND SIMULATION STUDY OF PROPERTIES OF ESTIMATORS 

The program described in Section II. B. 4 for computing the 

conditional density function in the GLAR1 process was used 

in two fashions. First, it was tested by using computer 

generated data from a GLAR(l) process with known parameter 

values k, q, and y . Simultaneously a simulation of properties 

^ A A 

of m. I.e ' s k and q for k and q was conducted. Second, it was 
used to estimate the parameters in the GLAR(l) model for real 
wind speed data. Only the first use is covered in this section. 
The second use is addressed in Chapter IV, Preliminary Data 
Analysis . 

Four aspects of the program were addressed in the use of 

the program with simulated data. These were: sensitivity 

of the maximum seeking method to start point, the size of the 

*maximum likelihood estimate 
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standard deviation of the m.l.e and moment estimates of k and 
q produced, and the degree of bias, if any, in the estimates. 

In addition, normality of the distributions of the estimates 
was investigated using normal plots and Shapiro-Wilks tests. 

Because of the large computation time involved in obtain- 
ing a m.l.e. the simulation study was small. Three types of 
data generated from GLAR(l) processes were used to exercise the 
program. Each type of data consisted of ten independent sets 
(replications) of 1000 data points each. The k and q values and 
the correlation were varied from one type of data to another. 
Thus the first type of data was generated with a k value of 
4.0 and a a value of 1.0. These parameter values produce a 
correlation of 0.75 (see equation II. B. 2.1). The second type 
of data varied the correlation, but retained the same k value. 

A k of 4.0 and a q of 3.0 produce a correlation of 0.25. These 
values were used to produce data sets of the second type. 

The third type of data returned to a high correlation, but 
used a small k value. The parameter values used to generate 
this data type were a k of 0.75 and a q of 0.1875. These 
values also produce a correlation of 0.75. Table II.E.l 
summarizes these cases. In all cases, y = 1. 

TABLE II.E.l 



CASE 


k 


q 


P 


1 


4.0 


1.0 


0.75 


2 


4.0 


3.0 


0.25 


3 


0.75 


0.1875 


0.75 



*maximum likelihood estimates 
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The Gamma variates were generated by the program LLRANDOMII 
[Ref. 16] and all runs were performed on the NPS IBM/3033 
computer . 

To test the sensitivity of the maximum likelihood computer 

program to the starting point of the search for a maximum, 

each set of data was used in two runs of the program. The 

first run used the actual parameter values k and q as a start 

point. The second run used the moment estimates k and q of 

the parameter values (see equations II. B. 4. 15 and II. B. 4. 16) as 

★ ^ 

a start point. The resulting m.l.e. parameter estimates k and 

/V 

q were recorded. 

The first case had k = 4.0 and q = 1.0. The results of 
the computer runs are presented in Table II. E. 2 for data of 
the first type. The last column presents the two-dimensional 
distance in the (k,q) plane between the estimates produced by 
the two different start points for each data set. Although 
the values do not differ widely, the relatively large differ- 
ences in some cases indicate that the likelihood function is 
relatively flat near the maximum. 

However, there is another factor which may be contributing 
to differences in final parameter estimates. The calculation 
of the likelihood function for 1000 data points requires the 
numerical evaluation of 999 integrals (see equations II.B.4.2 
and II. 3. 4. 3). The IMSL routine D CAD RE was used for this 
evaluation. Two of the parameters in the call to DCADRE are 



ie 

maximum likelihood estimates. 
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TABLE II. E. 2 



Results of Search for Maximum Likelihood Estimates 
in a GLAR(l) Model; k = 4.0, q = 1.0 



Data 

Set 


Starting 

k 


Value 

q 


Run Time 
(in Min.) 


Nunber of 
Iterations 


Ending Value Maximum 
k a Likelihood 


0 


4.000 


1.000 


20 


6 


3.668 


0.882 


-214.479 


0 


2.884 


0.538 


129 


17 


3.624 


0.886 


-214.488 


1 


4.000 


1.000 


19 


2 


4.004 


0.979 


-204.819 


1 


4.235 


1.101 


46 


5 


4.014 


0.983 


-204.820 


2 


4.000 


1.000 


104 


14 


3.451 


0.869 


-260.332 


2 


3.360 


0.805 


112 


14 


3.440 


0.867 


-260.333 


3 


4.000 


1.000 


16 


2 


3.977 


1.051 


-206.416 


3 


3.767 


0.947 


71 


15 


3.955 


1.041 


-206.418 


4 


4.000 


1.000 


52 


5 


4.246 


1.232 


-304.365 


4 


4.423 


1.233 


45 


6 


4.254 


1.235 


-304.366 


5 


4.000 


1.000 


47 


3 


4.061 


0.998 


-211.074 


5 


4.647 


1.187 


27 


4 


4.122 


1.028 


-211.060 


6 


4.000 


1.000 


61 


4 


3.593 


0.891 


-229.087 


6 


3.387 


0.782 


55 


9 


3.565 


0.883 


-229.092 


7 


4.000 


1.000 


37 


5 


4.041 


0.973 


-241.637 


7 


4.421 


1.069 


82 


7 


4.051 


0.974 


-241.636 


8 


4.000 


1.000 


37 


6 


4.024 


0.999 


-240.432 


8 


3.927 


0.856 


24 


8 


4.106 


1.033 


-240.424 


9 


4.000 


1.000 


51 


3 


4.109 


1.054 


-255.245 


9 


4.398 


1.112 


38 


5 


4.174 


1.081 


-255.229 



Value of 
Difference 
(10~ 3 ) 

46 



10 



11 



24 



8 



68 



29 



10 



88 

70 



88 



the relative and absolute errors allowed for this calculation. 

Practical considerations of computer run time dictated rela- 

-4 

tively large values of 10 for each of these parameters. 

This error when applied 999 times in the calculation of the 
likelihood function may have lead to the differences in m.l.e. 
parameter estimates. As a test of this hypothesis the data 
set which produced the largest distance between the pairs of 
estimates (data set 8) was rerun with DCADRE error parameters 
set at 10 The run which used the actual parameter values 

as a start point ran in 171 minutes and produced estimates of: 

A A 

k = 4.102, q = 1.031. The run which used the moment estimates 
as a start point was terminated after 404 minutes CPU time. At 

/\ A 

that point it had estimates of: k = 4.088, q = 1.025. The 

distance of 15 x 10 ^ represents a significant reduction in 

-3 

the previous distance of 88 x 10 . It seems from this exam- 
ple that the program can be made less sensitive to the starting 
point by increasing the accuracy with which DCADRE computes 
the integrals in the likelihood function. Of course, a con- 
siderable increase in computational cost is incurred. This 
cost is not practical in these simulations or necessary since 
only a rough idea of the properties of the estimates was sought. 
The results of the runs for data type one are also presented 

^ A A 

in Figure II.E.l. First, each pair of estimates, (k,q) and (k,q) 
is plotted in the k,q plane. Then each point is projected 
along each axis to more conveniently reflect the marginal 
variation in the estimates for each parameter. 
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E. (+) AND MOMENT ESTIMATES (CD) 

TRUE K = 4. 0 ; TRUE Q=1.0 
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When k = 4.0 and q = 1.0, the method of moments produces 
estimates for k with a sample mean of 3.94 and a sample 
standard deviation of 0.55. The statistics for the estimates 
of q are a sample mean of 0.96 and a sample standard deviation 
of 0.21. The maximum likelihood method produces estimates of 
k with a sample mean of 3.93 and a sample standard deviation 
of 0.27. The values for q estimates are a sample mean of 1.00 
and a sample standard deviation of 0.11. Although the method 
of moments and maximum likelihood method do not produce signi- 
ficantly different values for the mean of the estimates of the 
parameters, the lower estimated standard deviation for the 
likelihood method makes this technique more desirable from 
the standpoint of precision. No bias was evident in either 
estimation technique with the precision attained in the 
simulations . 

The second case was the low correlation case with k = 4.0 
and q = 3.0. Here the distinction between the two estimation 
procedures is not as clear (Table II. E. 3 and Figure II. E. 2). 

The method of moments produced estimates for k with a sample 
mean of 3.99 and sample standard deviation of 0.17. The esti- 
mates for q had a sample mean of 2.98 and sample standard 
deviation of 0.18. The maximum likelihood method produced 
estimates of k which had sample mean 4.00 and sample standard 
deviation 0.16. The q estimates had a sample mean of 2.99 
and sample standard deviation of 0.17. It is clear that neither 
method of parameter estimation enjoys a distinct advantage 
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TABLE II. E. 3 



Results of Search for Maximum Likelihood Estimates 
in a GLAR(l) Model; k = 4.0, q = 3.0 



Value of 



Data 


Starting 


Value 


Run Time 


Number of 


Ending 


Value 


Maximum 


Difference 


Set 


k 


q 


(in Min.) 


Iterations 


k 


q 


Likelihood 


(10~ 3 ) 


0 


4.000 


3.000 


38 


6 


4.078 


3.187 


-620.637 


8 


0 


4.019 


3.271 


25 


5 


4.085 


3.192 


-620.637 


1 


4.000 


3.000 


23 


4 


3.780 


2.934 


-655.488 


5 


1 


3.718 


2.843 


26 


4 


3.775 


2.931 


-655.489 




2 


4.000 


3.000 


40 


3 


3.789 


2.765 


-663.886 


3 


2 


3.912 


2.862 


27 


5 


3.789 


2.762 


-663.886 




3 


4.000 


3.000 


44 


7 


4.321 


3.132 


-585.776 


5 


3 


4.000 


3.116 


50 


4 


4.316 


3.130 


-585.774 




4 


4.000 


3.000 


44 


4 


4.166 


3.403 


-604.719 


22 


4 


3.958 


3.174 


68 


6 


4.152 


3.385 


-600.715 


5 


4.000 


3.000 


27 


6 


3.890 


2.904 


-638.667 




5 


4.001 


3.048 


13 


4 


3.885 


2.902 


-638.666 


5 


6 


4.000 


3.000 


43 


3 


4.051 


2.873 


-599.581 




6 


4.221 


2.902 


49 


5 


4.049 


2.878 


-599.581 


5 


7 


4.000 


3.000 


23 


2 


3.954 


2.934 


-627.062 




7 


4.079 


3.074 


16 


7 


3.963 


2.941 


-627.058 


11 


8 


4.000 


3.000 


21 


4 


3.917 


2.898 


-637.234 


11 


8 


3.708 


2.633 


32 


6 


3.908 


2.891 


-637.233 


9 


4.000 


3.000 


15 


2 


4.111 


2.924 


-590.563 


14 


9 


4.093 


2.906 


43 


2 


4.121 


2.934 


-590.562 
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over the other with respect to precision or accuracy. How- 
ever, the method of moments is considerably cheaper with 
respect to computation time required. This is consistent 
with known results that for i.i.d. Gamma data (k = q) , the 
moment estimate for k is quite efficient when compared to 
the maximum likelihood estimator. 

Third case (k = 0.75, q = 0.1825) . The third type of data 
was a high correlation case with a low k value. Specifically 
k = 0.75, q = 0.1825, and the correlation was 0.75. As can 
be seen in Figure II. E. 3 and Table II. E. 4, both the methods of 
parameter estimation considerably overestimated the parameter 
values, indicating considerable bias in the procedures. The 
method of moments produced estimates for k with sample mean 
of 0.8061 and sample standard deviation of 0.067. Those for 
q had a sample mean of 0.232 and sample standard deviation of 
0.026. The likelihood method produced estimates for k with a 
sample mean of 0.852 and sample standard deviation 0.039. The 
corresponding statistics for q estimates are 0.224 and 0.004. 
As was true in the other high correlation case, the standard 
deviations of the maximum likelihood estimates are consider- 
ably smaller than those of the moment estimates. However, 
since the evidence is that the estimates are highly biased, 
the advantage of this smaller standard deviation is not clear 
unless additional data would serve to reduce the apparent 
bias. The detailed results for this data type are presented 
in Table II. E. 4 and Figure II. E. 3. It would be of interest 
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TABLE II. E. 4 



Results of Search for Maximum Likelihood Estimates 
in a GLAR(l) Model; k = 0.75, q = 0.1875 



Value of 



Data 


Starting 


Value 


Run Time 


Number of 


Ending 


Value 


Maximum 


Difference 


Set 


k 


q 


(in Min.) 


Iterations 


k 


q 


Likelihood 


(10-3) 


0 


0.7500 


0.1825 


85 


5 


0.861 


0.221 


325.685 




0 


0.9795 


0.2489 


136 


4 


0.863 


0.222 


325.685 




1 


0.7500 


0.1825 


59 


4 


0.832 


0.235 


866.216 


n 


1 


0.8026 


0.2351 


19 


1 


0.832 


0.235 


866.214 


u 


2 


0.7500 


0.1825 


72 


4 


0.806 


0.224 


439.198 


1 


2 


0.7621 


0.2147 


66 


4 


0.807 


0.224 


439.198 


_L 


3 


0.7500 


0.1825 


141 


2 


0.866 


0.222 


1131.684 




3 


0.7514 


0.2213 


121 


1 


0.866 


0.222 


1131.671 


0 


4 


0.7500 


0.1825 


59 


4 


0.807 


0.223 


2333.130 




4 


0.7749 


0.2420 


54 


5 


0.807 


0.222 


2333.131 


1 


5 


0.7500 


0.1825 


59 


4 


0.910 


0.217 


883.801 




5 


0.8336 


0.2320 


210 


4 


0.910 


0.217 


883.801 


0 


6 


0.7500 


0.1825 


129 


5 


0.885 


0.222 


287.798 




6 


0.7308 


0.1790 


88 


4 


0.884 


0.222 


287.799 


1 


7 


0.2500 


0.1825 


71 


5 


0.795 


0.223 


977.265 




7 


0.7759 


0.2241 


117 


8 


0.795 


0.223 


977.265 


0 


8 


0.7500 


0.1825 


66 


6 


0.906 


0.227 


1018.470 


0 


8 


0.8092 


0.2863 


136 


6 


0.906 


0.227 


1018.472 


9 


0.7500 


0.1825 


57 


4 


0.818 


0.212 


410.842 


36 


9 


0.8335 


0.2347 


85 


3 


0.852 


0.225 


411.703 
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to see if a technique such as (2-fold) jacknifing, such as 
that applied by Quenouille [Ref. 17] to correlation estimates, 
would help here. 

A much larger study would be needed to come to definite 
conclusions about the efficiency of maximum likelihood esti- 
mation in this model. However, as in the i.i.d. case for 
Gamma variates, m.l.e. 's are likely to be considerably better 
than moment estimations for values of k less than one. 

Note too that the maximum likelihood estimation did not 
include the mean value parameter. This could be done or the 
mean could be estimated from the sample mean X. The inflation 
of variation of X due to the correlation is known to be 
(asymptotically) 



1 + 2 



CO 



l 

k=l 



k 

P 



1 + P 
1 - P 



Thus for p = 0.75, the effective sample size for estimating 
U in a sample size of size n is n (1-p) / (1+p) . For p = 0.75 
this is n/7. 

The normality of the estimates was investigated with normal 
plots and Shapiro-Wilk tests for normality. A summary of 
results if given in Table II. E. 5. The normality hypothesis 
is accepted at a 0.95 level if the Shapiro-Wilk statistic W 
is higher than 0.842, at a 0.99 level level if it is higher 
than 0.781. No strong indication of non-normality is indicated 
in any case. 



97 



Summary of Simulation Results 
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III. MOVING AVERAGE MODELS 



A. INTRODUCTION 

Although several researchers have proposed models for 
correlated, marginally exponential random sequences [Refs. 18, 
19,20], Gaver and Lewis [Ref. 2] produced the first analytically 
and computationally tractable model for the generation of 
correlated, marginally Exponential random sequences. They 
showed that in the usual linear, additive, first-order, auto- 
regressive equation 

X = 8X + E (III.A.l) 

n n+1 n 

where (X , n= 0,1,2,...} is a second-order stationary, 
n 

marginally Exponentially distributed sequence of random varia- 
bles; (E n , n = 1,2,...} is an independent, identically dis- 
tributed (iid) sequence of innovative random variables; 

0 < 8 < 1, the distribution of the (E } which produces the 
— n 

desired marginal distribution for {X n } is 

( 0 with probability 6 , 

E = ) (III. A. 2) 

n | E n with probability 1-8, 

where {E n , n = 1,2,...} is an iid sequence of Exponentially 
distributed random variables with the same parameter, A, as 
the {x^} sequence. Equation III.A.l can now be written as 
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with probability 3 , 



(III .A. 3) 



/ 3X , 

I n-1 

X = 
n 

' $X n _i + E n with probability 1-3. 

If (I , n = 1,2,...} is an iid binary sequence independent of 
(X n > and (E n > such that P(I n =0) = l-P(I n =l) = 3, then 
equation III. A. 3 can be more succinctly written as 



X n = 3X n-l + W (III.A.4) 

X n is a random linear combination of identically distributed 
random variables in the sense that the variable E n actually 
enters into the sum only when the random variable I has value 
one. Since for a given {i^} sequence, the distribution of X^ 
depends only on the distribution of X . and E , X will be 
Exponentially distributed whenever both X^_^ and E n are inde- 
pendent and have the same Exponential distribution. This 
understanding allows the autoregressive relation in III.A.4 
to be transformed into a moving average by substituting another 
innovative random variable for the to produce 

x n “ SE n + I n E n+1 . (III. A. 5) 

This model was designated the EMA(l) for Exponential moving 
average of order one. This EMA(l) model is one dependent in 
that X n and X n+ ^ are independent for j ^ ±1. Consequently, 
only the lag one correlation, = p_^ = CORR(X n ,X n+ ^) , or more 

completely only the joint distribution of X n and X^ + ^ need be 
studied. 
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Lawrance and Lewis [Ref. 5] give a fairly complete des- 
cription of this EMA(l) model. Of particular note is the 
relative tractability of this model which enabled the authors 
to derive correlations, distributions of sums of X^'s, 
intensity function, spectrum of counts, joint density of 
X n and X n+ ^, conditional expectations, and other properties. 

The existence of these characteristics is beneficial in data 
analysis and is a primary advantage of the EMA(l) over pre- 
viously suggested models. However, the EMA(l) model does 
possess a limitation. The range of possible positive correla- 
tions, p^, is restricted to the interval from zero to one- 
quarter. Thus for a given correlation between zero and 
one-quarter, the structure of and X n+ ^ and the sample path 
behavior of the sequence are determined. 

The structure of III.A.l is that of a special random linear 
combination of Exponential random variables to given an Exponen- 
tial random variable. Other such random linear combinations 
are now known and for the first-order autoregressive case 
produce dramatic differences in sample path behavior of the 
sequence (X n ) . In this section of the thesis we investigate 
these random linear combinations in the context of the first- 
order moving average structure. 

In this Chapter we examine extensions of the model in four 
ways : 

1 . Negative Correlation 

McKenzie [Ref. 21] has suggested a scheme for inducing 
negative correlation in the EMA(l) process by correlating the 
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sequence {I } . A better scheme, in that it produces a larger 
range of negative correlations, was introduced by Lawrance 
and Lewis [Ref. 8] for the extended first-order autoregressive 
model NEAR(l). This scheme involves a bivariate error se- 
quence {E n ,E^} and its use is investigated in this thesis for 
the moving average process . 

2. New Exponential Moving Average Model of Order One, 

NEMA ( 1 ) 

It will be shown that no first-order moving average 
process which is a random linear combination of Exponential 
random variables can have correlation greater than one-quarter. 
Thus the differences in the processes for given correlation is 
investigated in terms of the joint structure of X^ and x n+ ]_ • 

The NEMA(l) process obtained by using the NEAR(l) structure 
[Ref. 8] in a moving average context is analytically tractable 
and quantities such as the joint Laplace-Stilt jes transform 
of X n and x n+ j_/ the spectrum of counts, PCX^^ > , and condi- 

tional expectations are obtained. It also combines the forward 
and backward EMA(l) models as extreme cases and is thus a 
natural extension of the EMA(l) model. 

3 . The Moving Minimum Model 

A non-linear combination of Exponentials involving 
minima has been applied by Tavares [Ref. 22] to obtain a first- 
order autoregressive process which is intimately connected 
[Ref . 23] with the EAR(l) process of Gaver and Lewis [Ref. 24]. 
This structure is applied to the moving average process. It 
is shown that this process extends the range of attainable 
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correlations in first-order "moving average" process beyond 
that obtainable by random linear combinations of Exponentials. 
However, the process is slightly degenerate in that there is 
a positive probability that successive points will lie on the 
line X n+1 = bX n for positive correlations and on the curve 
e + e =1 for negative correlations. The price paid 

for the extended range of correlations is a limited analytical 
tractability as compared to the EMA(l) or NEMA(l) processes. 

The moving minimum process is investigated in terms of the 
joint structure of X n and X n+ ^ . Although the joint distribution 
can be derived, the functions are difficult to examine in 
detail. Therefore, simple characterizations of the joint 
structure, in addition to the correlation, are examined. These 
include the P(X n+ ^ > X^) , a crude measure of the tendency of 
the sequence to exhibit runs up and down, and conditional 
expectations, E(X n |X n _ 1 = x) and E (X n _^ [ X R = y) . 

4 . The Beta-Exponential Model 

Finally, another random linear combination of Exponen- 
tials to produce correlated Exponentials is examined. Unlike 
the previous models, the coefficients of the Exponential random 
variables are themselves continuous random variables. This in- 
creases the complexity and reduces the analytical tractability 
of this model. Simple sample path characteristics are derived 
or simulated. These are special cases of the GIMA(l) from Chapter II. 

B. NEGATIVE CORRELATION IN MOVING AVERAGE MODELS 

The problem of non-negative correlations was addressed by 
McKinzie [Ref. 21] who modified the form of the EMA(l) model to be: 
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(Ill .8.1) 



X 



n 




where (E n , n = 0,1,2,...} is an iid sequence of Exponential 
random variables, {Z , n = 1,2,...} is a sequence of binary 
random variables with P(Z n =l) = 1 - P(Z n =0) - 6, and Z n 
is independent of all E n and all past X^. By imposing an 
MA(1) correlation on the sequence {Z n }, McKenzie was able to 
produce a negative correlation for the {X^}. However, this 
negative correlation is achieved at the cost of reducing the 
possible range of positive correlations for the {X^}. Using 
McKenzie's formulation, the range of correlations obtainable 
with the EMA(l) model is (--gjr jgO • 

An alternative procedure for producing negative correla- 
tions was introduced by Lawrance and Lewis [Ref. g] . Their 
procedure requires two series of innovative factors {E , 
n = 0,1,...} and (E^, n = 0,1,...} which are correlated and 
may be antithetic. 

Antithetic variables are generated by using the fact that 
the variables E. and E.' defined as: 



where (U^, i = 1,2,...} is a sequence of uniform (0,1) random 
variables, are both marginally Exponentially distributed and 
correlated. P.A.P. Moran [Ref. 25 ] has determined that the 



l 



l 




-ln(U i ) 



(III .B.2) 



E[ = -ln(l-U i ) 
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correlation between E^ and is approximately -0.6449 

and this is the lowest correlation possible between 

Exponential random variables. E. and E! have a deterministic 

i i 

relationship since 



E! = 

l 



-E. 

-In (1 - e x ) . 



(Ill .B.3) 



Using the Lawrance and Lewis extension of the EMA(l) 
process, the model becomes 



= 3E + IE' . , (III.B.4) 

n n n n-1 

where {E n , n = 1,2,...} is an iid sequence of Exponential 

random variables, (E', n = 0,1,...} is a sequence of Exponen- 

n 

tial random variables which are correlated with the respective 
variables in the {E^} sequence, {I , n = 1,2,...} is a sequence 
of independent binary random variables with P(I n =0) = 1 - P(I n =l) 
= 3, 0 <_ 3 <_ 1, and {l n }, {E n } are independent of each other 
and all previous X n values . 

The correlation of the X's can then be computed as follows. 

Let E(X) = u and recall that since {X } and (E } have the 

n n 

same marginal exponential distributions, VAR(X^) = VAR(E n ) . 



X n+l X n - [BE n+l + 



= 3 E ,,E +61 ,,E E'+3I E ,.E' ,+I , i I E'E' , 

n+1 n n+lnn n n+1 n-1 n+1 n n n-1 
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Thus, using the independence of {E n >, { I n } , and the iid nature 

of {E } and {I } 
n n 

E (x n+ i x n ) = 6 2 y 2 +6(l-8) [COV(E n ,E^)+y 2 ]+6(l-6)y 2 +(l-6) 2 y 2 

= y 2 + 8 (l-6)C0V(E n ,E^) 



Therefore , 



cov(x n+1# x n ) = 6 (1-8 ) COV (E r ,E^) 

and using VAR(X) = VAR(E) 

CORR(X n+1 ,X n ) = 8 (l-8)C0RR(E n ,E^) (III.B.5) 

Using Moran's result [Ref. 25] the correlation of antithetic 
Exponentials is known to be approximately -0.6449. Therefore, 
the greatest negative value that can be achieved for CORR(X n+ ^ ,X n ) 
is approximately -0.1612 when 8 = 0.5. Since no restriction 
was placed on CORR (E /E^) , the sequences (E n ) and {EJ^} need 
not be antithetic, but can have any correlation that is possi- 
ble for two Exponential sequences with the same marginal dis- 
tribution. By specifying that E^ = E , the original EMA(l) 
is achieved and the correlation for the X's is 8(1-6) as in 
Lawrance and Lewis [Ref. 5]. By allowing the correlation be- 
tween the ( E n } and {E^} sequences to vary from -0.6449 to 
1.0, the correlation of the X's will vary from a minimum of 
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-0.1612 to a maximum of 0.25 (also depending on the value of 
6) as can be seen from III.B.5. The Lawrance and Lewis exten- 
sion of the EMA(l) model gives greater possible variation in 
the correlation than that of McKenzie, but requires two se- 
quences of Exponential random variables to achieve this range. 

Although it is clear that with E' = E , a C0RR(E ,E') =1 

n n n n 

and when {E } and { E * } are antithetic CORR(E ,E') = -0.6449 , 
n n n n 

it may not be obvious how to generate (E^} sequences with 
correlations between these two extreme values. A simple bi- 
variate exponential sequence with any desired correlation in 
the permissible range can be generated using the relationship 

| E^ with probability p, 

E^ = (III.B.6) 

' E^ with probability 1-p, 

where E^ is the antithetic of E^. Then the extended EMA(l) 
model has two parameters, 6 and p, and the range of correla- 
tions for the X's is -0.1612 to 0.25, as above. The bivariate 
density for the pair {E^,E^} is not smooth. Other bivariate 
densities such as those in Gaver [Ref. 26] and Lawrance and 
Lewis [Ref. 27] can also be used. 

The above ideas on extending the correlation structure of 
the moving average models to encompass negative correlation 
can be applied to all of the new models given below. Details 
are not given. 
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C. THE EXTENDED EMA(l) MODEL , NEMA(l) 

1 . Introduction 

The original Exponential moving average process is 
discussed by Lawrance and Lewis [Ref. 5] . This paper 
considers the first order process defined by: 



X 



n 



SE n with probability 8, 

(III.C.l) 

SE n + E n+ ^ with probability (1-3), 



where {E n , n = 0,±1,±2,...} is an iid Exponential sequence and 

0 <_ 8 <_ 1. This random linear combination of Exponential 

variates is called the EMA(l) model for Exponential moving 

average of order one. Since X is a function of E and E , . , 

^ n n n+1 

this version is called the forward EMA(l) . The backward 
version of EMA(l) is obtained when E n+ ^ is replaced by E n-1 
in III.C.l. 

The fact that EMA(l) is a single parameter model sug- 
gests that this model may not be sufficiently flexible to 
address all processes. Investigation of an alternate, two 
parameter model may indicate that a two parameter model is 
sufficiently more flexible to justify its increased complexity. 

The extended, two parameter model is based on the new 
Exponential autoregressive process of order one (NEAR(l)). 

The NEAR(l) model propounded by Lawrance and Lewis [Ref. 8] 
is defined as 



103 



with probability a 



(III. C. 1.2) 



X 



n 




+ BX 



n-1 



with probability (1-a) 



where {X n , n = 1,2,...} is a second-order stationary sequence 
of Exponential random variables with parameter X, {E^} is an 
iid sequence of innovative factors, 0 _< B £ 1 , 0 <_ a <_ 1 , and 
a+B <2. By letting ( s ) and <p (s) denote the Laplace- 
Stieltjis transform of X and E respectively, Lawrance and Lewis 
determined that <j> £ (s) = ’ x+ (l-a) Bs * Us i n 9 a partial 

fraction solution technique to invert this transform produced 



1“ s 

! E n with probability (i- a ) g 

(III .C. 1.3) 

ctB 

( 1 -a) BE with probability ^ — jV- — 7-3-. 
n 1- ( l-a) B 

where {E n , n = 0,1,2,...} is an iid sequence of Exponentially 

distributed random variables which has the same parameter as 

the (X } sequence, 
n 

By noticing that the autoregressive model given in 
III. C. 1.2 using III. C. 1.3 is a random sum of two iid Exponen- 
tially distributed random variables, the NEMA(l) model is 
produced by substituting E n _^ for X n _^ in the NEAR(l) model. 

This procedure is identical to that used to produce the EMA(l) 
model from the EAR(l) model and yields 

( E n + SE n _^ with probability ct, 

X = (III. C. 1.4) 

n | 

' E^ with probability 1-a. 
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This model can be written in a more compact form as 



X 



n 




(III. C. 1.5) 



where (X^, n = 1,2,...} is a second order stationary sequence 
of marginally Exponentially distributed random variables; 

(E n , n = 0,1,...} is an iid sequence of Exponential random 
variables with the same parameter as the {X n } sequence; 

{I , n = 1,2,...} is an iid sequence of random variables with 
P(I = 6 ) = l-P(I n =0) = a; (K n , n = 1,2,...} is an iid se- 
quence of random variables with P (K =1) = 1-P(K = (l-a) 3) = 

n n 

I - 8 

]__ ' ( l-a ' ) 8 ; ^n^ 7 ^ K n^ 7 an< ^ ^ E n^ are rautua Hy independent; 
0<a<l;0<8<l; and a+8 < 2. 



ward versions of EMA(l) as special cases. When a = 1; 
P(I =6) = 1/ (l-a) 3 = 0, and P(K n =0) = 6. Hence, the 

NEMA(l) model becomes 



The NEMA(l) model contains both the forward and back- 



I 

I 




with probability 6 



X 



n 



(III .C.l .6) 



3E . + E with probability (1-8) . 
n-i n 



This is a form of the forward EMA(l) . 



When 8 = 1; P (I = 1) 



a, (l-a) 8 = (l-a), and 



P(K n = (l-a)) 




1. In this case, the NEMA(l) model 



becomes 
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with probability (1-a) , 



( 1-a) E 



n 



X 



n 



(III. C. 1.7) 



(l-a)E n + E n _^ with probability a 



This is a form of the backward EMA(l) with 8 = 1-a. There- 
fore, the NEMA(l) contains the special cases of EMA(l) when 
a and 8 assume specific values. 

One should also note in passing that the {X n > sequence 
becomes an iid sequence if a = 0 or 8 = 0. 

2 . Correlation Structure 

The following relationships will prove of value in 
succeeding calculations 




P(I n =0) = a. 



(III. C. 2.1) 




a8 + (1-a) *0 



a8 • 



(III. C. 2. 2) 



P(K n =l) 



l-P(K n = (1-a) 8) 



(III. C. 2. 3) 



1-8 



1- ( 1-a) 8 





) + ( 1-a) 6 ( 



1- (1-a) 3 



2 2 2 

1-3+ae -a B 

l-8+a8 



l-8+a8-a3+a8 2 -a 2 3 2 
1- 8+a6 



l~8+a8 _ a8 ( l-8+a6 ) 
l-8+a8 l-8+a8 
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E( V 



1-aB 



(III. C. 2. 4) 



X 



n 




(III. C. 2. 5) 



E(X ) 
n 




(III .C. 2 .6) 




E (K )E(E ) + E ( I ) E (E . 
n n n n-1 



(l-a3)E(E ) + a3E (E . 

n n-i 



E ( X) 



E (E) 



Since {X } and {E } are both Exponential, E(X) = E (E) implies 
n n 

VAR(X) = VAR(E) . Since both (E } and {X } have the same 

n n 

Exponential parameter, without loss of generality this param- 
eter will be considered to be 1. This, of course, requires 
E ( X) = 1 and VAR ( X ) = 1. 

The possible range of correlations for the NEMA(l) 
model can be determined by a simple calculation. We have 



X 




n 



X 



n+1 




Thus 
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X X . = (K E + I E .) (K , , E -+I , . E ) 

n n+1 n n n n-1 n+1 n+1 n+1 n 



= K , K E , , E + K ,.I E ^ , E . + K I E 
n+1 n n+1 n n+1 n n+1 n-1 n n+1 n 



+ I , , I E E , 
n+1 n n n-1 



Using the independence of {K }, {i } and (E } and the iid 
3 n n n 

nature of these three sequences, we have 



E(X n x n +i ) = (1-aS) 2 [E(E n ) ] 2 + ( 1-aS ) a6 [E (E^) ] 2 

+ (l-o8)a6 [2 {E (E n ) } 2 ] + (a6) 2 [E(E n ) ] 2 



= 1 + (l-aB)aB 



Therefore , 



COV(X n ,X n+1 ) = (l-aB)aB 

and 

CORR (X_ , X ) = (l-aB)aB. (III.C.2.7) 

n n+1 

Therefore, the original NEMA(l) model has the same range 
of possible correlations as the EMA(l) , namely that the corre- 
lations must lie in the interval [0,^-]. One should note that it 
is not possible to distinguish the parameters a and B from the 
correlation, or even whether the product, aBr has a given value 
or one minus that value. This is similar to the normal moving 
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average model of order one where the cases <p and 1 /$ are 

indistinguishable. In the normal case, the range of <j> is 

limited to the interval [0,1] on the basis of invertibili ty 

[Ref. 28] . It would seem simple and convenient here to limit 

a8 to the interval [0,-^]. However, non-normal processes are 

not completely determined by their correlation structure. In 

fact, Jacobs and Lewis [Ref. 6 ] showed that in the EMA(l) 

case, the values 8 and (1-6) can be distinguished using direc- 

2 2 

tional moments, E ( x n x n+ ]_) and E ( x n x n+ ]_)* Hence, such a re- 
striction on the value of ot6 is inappropriate. 

One should also note that, since the correlation 

between X and X , T , is zero for all K with absolute value 
n n+K 

greater than one, the first-order correlation completely 
determines the correlation structure of the model. 

The restriction on the range of attainable correlation 
is disappointing but not surprising since it can be proven that 
any Exponential moving average process of order one generated 
as a linear combination of independent Exponentials must have 
a correlation that lies in the range [0,-i], The proof of 
this contention follows the previous calculation closely. 

THEOREM: 

Assume {E n , n = 0,1,2,...} is a sequence of iid Exponen- 
tial variables with unit mean, {A n , n = 1,2,...} and {B r , 
n = 1,2,...} are sequences of iid random variables, and {A^}, 
{B^} , {E^ } are all mutually independent. Moreover, assume the 

sequence {X n , n = 1,2,...} defined by 
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(Ill .C. 2 . 8) 



X = A E + B E , 
n n n n n-1 



is a unit mean Exponential sequence. Now 



E (X ) = E ( A E + BE . ) 

n n n n n-1 



- E (A ) E (E ) + E(B) E(E n . ) 
n n n n- 1 



1 = E (A ) + E (B ) 

n n 



(III .C. 2 .9) 



In addition, since X n > 0 for all n, both A n and B p must be 

non-negative for all n. Hence E(A ) >0 and E(B ) > 0. It 

n — n — 

also follows that 1 > E (A) and 1 > E(B) . Now 



X X ^ 
n n+1 



= (A E +B E .) (A . E J _ 1 + B . E ) 
n n n n-1 n+1 n+1 n+1 n 



Therefore 



E ( X X . n ) 
n n+1' 



E (A , , A E ,.E + A , , B E . . E . + A B . . E 

n+1 n n+1 n n+1 n n+1 n-1 n n+1 n 



+ B , . B E E . ) 
n+1 n n n-1 



[E ( A) ] 2 + E (A) E (B) + 2E (A) E (B ) + [E (B) ] 2 



= [E (A) + E(B) ] + E ( A) E (B ) 



Since E (A) + E (B) = 1 from III. C. 2. 9, E (A) + E (B) = 1 



E(X n X n+1 ) = 1 + E ( A) E (B) 
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Therefore, since E(X n ) is one, by assumption 



COV ( X , X 



n ' n+ 1 



E (A) E (B) 



E (A) [1 - E ( A) ] 



and 



CORK (X ,X. 



n' n+1 



E (A) [1 - E (A) ] 



(III .C. 2. 10) 



Since 0 < E(A) < 1, the correlation must lie in the interval 



reformating the model. We choose to do this first by using 
the method devised by Lawrance and Lewis [Ref. 8] and given 
in equation III.B.4. With variables and sequences defined as 
in equation III. C. 1.5 let {E n } and {E^} be correlated sequences 
and redefine the NEMA(l) as 



[0,i]. Q.E.D 



The possible range of correlations can be extended by 



X. 



n 




(III. C. 2. 11) 



Then it follows that 





n+l E n+l + I n+l E n ) 
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Thus 




(1-aB) 2 + (l-aB)aB + (l-aB)aB [COV[E ,E']+1] + a 2 B 2 

n n 



1 + (1-aB ) aBCOV (E ,E') 

n n 



and 



COV(X ,X. 



n' n+1 



(1-aB) aBCOV(E n ,E^) 



Fina ly 



CORR(X ,X. 



n' n+1 



( 1-aB ) aBCORR (E n , E^) . 



(IIIC .2.12) 



As described above, Moran [Ref. 25] has determined that 



the range of possible correlations for two Exponentials is 

(-0.6449,1.0). Thus when aB = 0.5, the possible correlations 

for X and X fall in the interval (-0.1612,0.25). This 
n n+1 

procedure extends the range of possible correlations at the 
cost of generating the additional {E^} sequence. 

McKenzie [Ref. 21] has suggested that the range of 
correlations could be extended by requiring that the { I n } 
sequence be correlated. Using this scheme, he was able to 
show that the correlations for the {X n } sequence lies in the 
interval (- '^j/yg-) ♦ Because of the requirement of the moving 
average process of order one to have zero correlation for lags 
greater than one, the correlation of the (l n ) sequence also 
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had to be of MA(1) type. It is this restriction that pro- 
duces such a narrow range of correlations. A logical and 
obvious extension to McKenzie's work is to require that both 

the (K } and {I } sequences in the NEMA(l) model have a MA(1) 
n n 

correlation structure. This can be combined with the corre- 
lated {E }, {E'} scheme of Lawrance and Lewis. If this com- 
n n 

bined case is carried out, the NEMA(l) model is as follows 



X 



n 



K E 
n n 



IE' . , 
n n-1 ' 



(III .C. 2. 13) 



where (K n , n = 1,2,...} is a sequence of random variables 

with an MA(1) correlation structure with P(K n =l) = 

1“ 3 

l-P(K n = (l-a)B) = (i-ct ' ) ~ g • n = 1/2,...} is a sequence 

of random variables with an MA(1) correlation structure with 

P(I =8) = 1-P(I = 0) = a; {E } and (E'} are correlated se- 

n n n n 

quences with marginal Exponential distributions with unit 
means; and {K n }, (I }, and {E n } are mutually independent. 

Now 



X X = (KE + I E ' . ) (K . . E , -, + I ,iE') 

n n+1 n n n n-1 n+1 n+1 n+1 n 



K ,K E , . E + K ..I E ... E 1 ,+K I , . E E' 

n+1 n n+1 n n+1 n n+1 n-1 n n+1 n n 



+ I . I E ' E ' 
n+1 n n n-1 



So 
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E (X X ) 
n n+1 



= E (K K ) + ( 1-aB ) a8+ ( 1-aB ) aBE (EE') +E ( I. . I ) 
n+i n n n n+1 n 

= C0V(K n+1 ,K n )+(l-a8) 2 +(l-aB)aB 

+(l-aB)aB* [COV(E ,E')+1]+C0V(I ^,1 )+(aB) 2 

n n n+l n 

= l+COV(K n+1 ,K n ) +COV ( I n+1 , I n ) + (1-aB) aBCOV(E n ,E^) 



Therefore 

COV(X n> X n+1 ) = COV(K n+1 ,K n )+COV(I n+1 ,I n )+(l-aB)a6COV(E n ,E^) 
and 

CORR(X n ,X n+1 ) = COV(K n+1 ,K n )+COV(I n+1 ,I n ) (III. C. 2. 14) 

+ ( 1-aB ) aBCOV (E ,E') 
n f n 

Although this scheme obviously extends the range of possible 
correlations, it does so at the expense of considerable com- 
plexity. Considering the limited range of correlations 
possible by imposing a correlation on the {l n } and 
sequences, the additional complexity may not be warranted. 

If in spite of the complexities involved, one decides to in- 
duce correlations in the coefficient sequences, the NEMA(l) 
because it has two such sequences will yield a slightly larger 
range than the EMA(l) model. 
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3. 




One of the possible advantages of a two parameter 



model is the capacity for modifying p ( x n+ ]_ >x n ) and, conse- 
quently, the sample path behavior of the process while main- 
taining a constant correlation. Since the correlation is a 
function of a$, one can vary the values of a and 6 while 
keeping the product constant. The p ( x n+ ]_ >x n ) can calcu- 
lated by addressing each of the sixteen possible combinations 
of K and I values for x n+ ]_ and X n , computing the probability 
for each combination, and weighting the probability associated 
with a given combination by the probability that the given 
combination occurs. A sample calculation is provided and 
complete results presented in Table III. C. 3.1. 



Example: We have 



X. 



n 




X 



n+1 




P(K n =l) 



1-P (K n = (1-a) S) 



1-S 



1- (1-a) 3 



P(I n =3) 



i-p(i n = o) 



a . 




6 and consider the case where 



I 



n 



3, K 



n 



1, I 



n+1 




1 
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Since the [l n > and [K n > sequences are both iid and 

independent of each other, the probability of this combination 

of parameter values is simply the product of the individual 

probabilities of occurrence. Hence the probability of 

2 2 

occurrence is a 6 . Then in this case 



P(X nH.l >X n> “ P(E n+l +6E „ >E n +6E n-l> 



■ P (E , > (l-B)E +6E , ) 

n+1 n n-1 



(III. C. 3.1) 



Now is independent of Y = (1-8) E n + ^ E n _^ . Therefore, the 
calculation required by equation III. C. 3.1 is straightforward 
once the density of Y is obtained. We have, with f £ (x) 



-x. 



n-1 



the p.d.f. of E n-1 (i.e. e ) 



y/6 

P ( [1-8 ] E +SE i < y) = / P ( [1-6 ] E +8x < y I E . = x) f (x) dx 

n n-l- q n n-l 



y/6 

/ P( [1-8 ]E„ < y-SxjE^ . = X) f_ ( x) dx 
J r, n — ^ 'n-l E , 

0 n-1 



y/s 



/ P(E niW! E n-l = x > E E < x)d * 

0 n-1 



y/s -iff _ x 

/ ( 1 - e ' ) e x dx 



y/6 y/6 - U-ff 

J eax-J e e 
0 0 



dx 



_-JL_ y/s - :: r 1 ~ 2S 1 
i-e-y/S-e 1-6 J y/3 e X 1-6 dx 

0 
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P( [l-8]E n +6E n _ 1 < y) 



y/R _ x f j- ] 

1 - e-y/8 - e 1 * 6 [£* ] ^ [£«]. ^ 



dx 



= 1 - e 



- .-y/3 _ e "l-8 r Izl 



y r i-2g 1 

/ -I _ 6 1 1-8 \ 

^1-28^ (1_e ) 



= i- e ~Y/Z . rJ^Lie 1-6 + rJ^_i e -y /B 
1 e 1-28 e l I^28 Je 



8 r, -y/8 n A 1-6 1-3 



P(Y<y) = - j^[l- e * /iJ ] + 



1-23 



_-JL_ 

(1-e I'*) 



Therefore, the density function of Y, fy(y) , is 

- ( x ~ ~ 2 ~ 3 ) (-|) e y// ^ + ( j ~ Yg ) ' 3 e -y// ^ 1-1 ^ , for 6 ^ -j. Gaver and 

Lewis [Ref. 2] gave the necessary and sufficient conditions 

-X.x ~' a 2 x 

for a mixed exponential of the form iT^X^e +ir 2 X 2 e , 

i 1» + 71 2 ~ 1' an< ^ A 1 < ^2 to a P ro P er density 

function. The condition is that 



7T 



1 



< 




(III .C. 3 .2) 



In this situation we address two separate cases. The first 

1 g 

case is when 0 < 8 < -j. In this case - j -~ 22 > 1 and the 
requirement is 



1-8 

T=7Z 



< 



[1 ■ ^ire } v (f) >] 1 £ U 




-1 



1-3 

1-28 
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And III. C. 3. 2 is satisfied. The second case is when ^ < 8 < 1 

ft 

In this case - y~2 g > ^ an< ^ t ^ ie r equirement is 



- i i 1 - «!> * < i =*>>] 



-1 



< u-i^r 1 



< - 



1-26 



And III. C. 3. 2 is satisfied. If 6 = j, the density of Y is a 
Gamma(2) density. Therefore, the p.d.f. of Y is a proper 
density. This result can be used to complete the calculation 
of p ( x n+ ]_ >x n )* Recall that equation IIIC.3.1 stated 



P (X > X ) 
n+1 n 



= P (E , , > (1-8) E +6E , ) 

n+1 n n-1 



P(E n+ i>Y) 



By conditioning on Y and using the p.d.f. for Y derived above 
this is 



P (X > X ) 
n+1 n 



/ p (E n+1 > y I y = y) f y (y) d y 



= f e ' y [- ( 8 ) (h 
> ® 1 l l-28 ; V 



v 

^6 



— ' 1 'i> e -y /6 + (^.) (^(e * “)dy 



,6+1, 

f 1 w 6 v f°° ,3+1. " Y( 3 } 
1-28 ( 6+l J * 0 ( 8 } 



dy 



„ - 

1 - 6 , f , 2 - 6 , ^ 1 - 6 ' 



i 1 , , 1 ~ 6 , ( ,2 — 8 , „ 

1-26 2-8 > 0 ( l -6 



dy 
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P(X n+l >X n> 



-( 



) <: 



1-23 7 V S+1 



-) + ( 



_i_) (1=8) 

1-23 ^ ^2-3 7 



P (X ... > X ) 
n+1 n 



(2-3) (1+3) 



(III. C. 3. 3) 



Table III. C. 3.1 presents the results of the calculations for 

all of the sixteen combinations of parameter values for X n+ ^ 

and X . 
n 

When the P (X > X ) was calculated for various values 
n+l n 

of a and 3/ it was found that the values for this probability 
varied from 0.44 to 0.54. Table III .C .3 . 2 contains the results 
of these calculations for four hundred forty-one combinations 
of parameter values. Although the variation in probability is 
not large, it does represent an increase over the forward 
EMA(l) model. In particular, the forward EMA(l) model can not 
produce a probability greater than 0.50. Consequently, the 
NEMA(l) model not only has a greater range of possible proba- 
bilities, but also can produce probabilities greater than 0.50. 
The implications of this greater range is that the NEMA(l) model 
can address data sample paths that have a slight tendency for 
either runs of increasing or decreasing values, while the 
EMA(l) can only address sample paths that tend to produce runs 
of decreasing values. 

Examples of scatter plots and sample paths for three 
sets of parameter values and positive correlations are given 
in Figures III .C . 3 . 1-III.C . 3 . 6 . Because of the relatively low 
correlations possible and because of the limited range of 
values for P ( x n+ j_ > X^) , differences among the figures are not 
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TABLE III. C. 3.1 
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sharply delineated. However, differences can be detected, 
particularly in the sample paths. In Figure III. C. 3.1 with 
an a value of 0.95 and 6 value of 0.50 has a P(X ,, > X ) = 
0.44, the lowest value for this probability. Since this pro- 
duces a slight tendency for runs of decreasing value, the 
number of extreme values (i.e. greater than 3.0) is two. 

In Figures III. C. 3. 2 and IIIC.3.3 the P(X > X ) is 0.50 and 

n+i n 

0.54, respectively, with a corresponding increase in 
the number of large values. This trend is more 
difficult to detect in the corresponding scatter plots. 
Figures III .C . 3 . 7-III . C . 3 . 12 provide sample paths and scatter 
plots for the same a and 8 values as previously displayed, 
but with antithetic innovative sequences (see III.B) and 
consequent negative correlations. Although the negative 
correlation is evident, trends in these figures are difficult 
to detect. The extremes of sample path variability produced 
by the NEAR(l) process [Ref. 8] are not reproducible with the 
NEMA(l) process. This may be attributable to the restricted 
range of possible correlations. 

4 . Laplace Transform of Sums 

One aspect of the EMA(l) model is its analytical 
tractability . This is evidenced by the ability to derive the 
Laplace transform of sums by a recursive relationship given 
by Lawrance and Lewis [Ref. 5] . This tractability carries 
over to the NEMA(l) process. The Laplace transform is useful 
in obtaining quantities which are of use in analyzing point 
processes, namely the intensity function and the spectrum of 
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FIGURE III. 
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FIGURE III. 
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FIGURE III. 
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FIGURE III. 



EM A SAMPLE PATH 
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FIGURE III. 
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counts. These quantities are derived for the NEMA(l) process 
in subsequent sections from the results obtained here. 



1-6 



With X defined as in III. C. 1 . 5 , let — -n — 6__ = 5 
n ' 1- (1-a) 6 



-sT 



and (1-a) 6 = y. Further, let £ X. = T and <J>_ (s) = E(e r ) . 

1=1 r 



Then we have 



T = X. + X_ + . . . + X 
r 1 2 r 



(III .C . 4 . 1 ) 



K^i + i iE q + ^ 2^2 + ^ 2^*1 + • • • + ^r^r + -*- r E r -i 



K E + (I +K . ) E . + ... + ( I- + K. ) E. + I.E n 

rr r r-1 r-1 2 11 10 



Then letting L ^ = Ij + ^ + Kj, j = l, 2 ,...,r-l and using the 
mutual independence of the iid sequences {K n >, {l n l. 



-sT r 

4> t (s) = E (e ) 



(III .C. 4 . 2 ) 



— s[K E +L .E . + . . . +L. E. + 1 . E-. ] 

_ , r r r-1 r-1 11 10, 

E (e ) 



-sK E -sL . E . -sL, E, -sI.E n 

= E(e r r )E(e r ‘ 1 ...E(e ^(Efe 10 ) 



-SK.E. -Sl.E. -SL.E. 

Now let 'F„(s) = E(e ^ , ^(s) = E(e ^ ^) , (s) = E(e -* 



K 



Then 



4> t (s) = T k (s)'F i (s) [Y l (s) ] 



r-1 



(III .C. 4 . 3 ) 



To evaluate these quantities note that 
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K = 



1 with probability 6 , 



y with probability 1-6 



Then 



-SK E. 

Y k (s) = E(e 3 3 ) 



-sE. -syE. 

= 5E(e - 1 ) + (1-6 ) E (e 3 ) 



V K (s) = 6<|» (s) + (1-6)<J> (ys) , 



(III .C. 4 .4) 



where <J> E (s) = 



So 



V S) 



6 + ( 1 - 6 ) 



1+s 1+ys 



(III. C. 4 .5) 



Also 



I = 



8 with probability a. 



0 with probability 1-a, 



so that 



-sI.E.) -s6E. 

'i' (s) = E (e 33 = aE (e 3 ) + (1-a) 



(III. C. 4 .6 



^1 ^ l+8s + ^ 1-a ^ • 



(III .C .4 . 7) 
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To evaluate ^(s) note that 



L = 



3+1 


with probability 


a6 , 


3+y 


with 


probability 


a (1-6 ) , 


1 


with 


probability 


(1-a) 6 , 


Y 


with 


probability 


(1-a) (1-6) . 



Therefore 



Y l (s) = E(e 3 3 

-S [3 + 1] E . -s[3+y]E. 

= a6E (e J) + a (1-6) E (e 3 ) 

-sE. -syE . 

+ (l-a)6E(e 3 ) + (1-a) (l-6)E(e 3 ) 



Y l (s) = a6<J> E ( [0+l]s) + a(l-6) <t> E ( [0+Yls) + (l-a)6<J> E (s) (III. C. 4. 8) 

+ (1-a) (1-6 ) <J> E (ys) . 

Using the results of III. C. 4. 3, III. C. 4. 4, III. C. 4. 7, and 
III. C. 4. 8 



<J> T (s) = [od E (s) + (1-6 ) 4> E (ys) ] x [atf E (3s) + (1-a) ] (III. C. 4. 9) 

r 

x [aocj> E ( [3+D s) +a (1-6) <j> E ( [3+y ] s) + (1-a) 64> E (s) 

+ (1-a) ( 1-6 ) 4> E (ys ) ] r_1 
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This extends the result (3.7) in Lawrance and Lewis (Ref. 5] 

to the NEMA(l) process. 

5 . Laplace Transform of the Distribution 
of ~ Counts 

The Laplace transform of the sum is useful in deriving 
the distribution of the synchronous counting process of the 
number of events that occur in ( 0 , t ] when the origin is estab- 
lished at the occurrence of an arbitrary event. The number 
of events in (0,t] is related to the distribution of a sum by 
the relationship 

n 5 < r iff T > t, r = 1,2,... (III. C. 5.1) 
t r 

where n£ is the number of events in (0,t] and T is the sum 
t r 

of the first interevent times. Thus 



P(N^= r) = F r (t) - F r + 1 (t) 



where F (t) is the distribution function of T . The proba- 
r r 

bility generating function of can then be written as 



(z;t) 




) 



CO 

= l Z r P(N^ = r) 
r=0 r 

co 

= l Z r [F r (t) - F r+1 (t)] 
r=0 

CO 

= 1 + (Z-l) l Z r_1 F (t) . 

r=l 
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* * 
Let 'F f (z;s) be the Laplace transform of H' f (z;t) r and f r (t) 

the Laplace transform of f r (t), the p.d.f. of T . Then 



Y*(z;s) 



1 

s 



(1-Z) 

s 



00 

r r*— 1 * 

1 z r 1 f r (t) 

r=l r 



(III. C. 5. 2) 



Using the Laplace transform of the sum from III. C. 4. 9 



Y * f ( z ; s ) 



V * ( z ; s ) 

X 



uu 

" - ( -~ 2) I z r_1 [ 0 (j) E (s) + (1-6 ) <P E (ys) ] 



r=l 



x[a<(» E (es) + (l-a) ]x [a6<^ E ( [0+l]s)+a (1-6 ) <|> E ( [3+Y ] s) 



+ (1-a) 6<j) £ (s) + (1-a) (1-6 ) <j> E (ys) ] 



r-1 



-^^[6<J> e (s) + (1-6)<D e (ys) ]x[a6 £ (3s) 



(III. C. 5. 3) 



+ (l-a) ] x< 



1-z [a6<t> E ( [3+1] s+a (1-6 ) <j> E ( [0+y] s) 
^+(l-a)66 E (s) + (l-a) (1-6 ) <|> E (ys) ] 



where <D E (s) = 

If m^(t) is the intensity function of the point proc- 
* 

ess, then m f (t), its Laplace transform, can be obtained by 
differentiating III. C. 5. 3 with respect to z, evaluating the 
derivative at z = 1 , and then differentiating with respect to 
s. These steps, when taken, produce a series of tedious 
calculations which produce no analytical insights. The result 
of these steps is 
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l+(l+46-2a6) s+( 36+56 2 -2a6-5a6 2 +a 2 6 2 ) s 2 



m. 



2 3 2 322233 

+ (26 +26 -3ctB -3ag J +ct g z +ct B ) s J 

s+(l+46-3a6+a 2 6 2 ) s 2 +(36+56 2 -2a6-7a6 2 +3a 2 6 2 ) s 3 



(III .C.5 .4) 



2 3 2 322214 

+(26 +28 -3a6 -3a8 +a 8 +a 8 )s 



This result can be verified in a number of ways. First, when 

a = 1, the process is the EMA(l) process and, hence, the formula 

must reduce to (4.2) given in Lawrance and Lewis [Ref. 5] with 

A = 1. Second, when a = 0, the NEMA(l) process reduces to a 

Poisson process and the formula under this condition must 

reduce to the Laplace transform of the constant intensity 

function of a Poisson process with rate 1, Third, with 

8=0 the NEMA(l) process is again a Poisson process. Finally, 

* 

using one of the Tauberian Theorems, lim m f (t) = lim sm, (s) = 1. 

t+<» 1 s->0 r 

We take these cases in turn. First, when a = 1 



l+(l+46-2a8) s+(36+58 2 -2a6-5a6 2 +a 2 8 2 )s 2 



* 

m f (s) 



2 3 2 322233 

+(28 +28 -3a6 -3aB +a B +a B ) s 

s+(l+46-3a6+a 2 6 2 )s 2 +(36+56 2 -2a6-7a6 2 +3a 2 6 2 ) s 3 



2 3 2 322234 

+(28 +28 -3a8 -3o8 +o 6 +a 8 ) s 



reduces to 



1+(1+2S) s+(B+B 2 ) s 2 
s+ (1+B+B 2 ) s 2 + (B+8 2 ) s 
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t 



m* ( s) = ( 1+6 s ) ( 1+ [1+6 3 s ) 

f s( [1+s] [1+(B+B 2 )s] ) 

= (l+6s) (1+ [1+6 ] s) 

6 (1+6) s (1+s) ( g~+gj + s) 

which is the result of Lawrance and Lewis [Ref. 5] with X = 1. 
In the second case with a = 0 

_ 1+(1+46)s+(36+56 2 )s 2 +(26 2 +26 3 )s 3 

m f ' S ' 2 2 3 2 3 4 

1 s+ (1+46) s + (33+56 ) s J + (26 +26 )s* 

1 

s ' 



the Laplace transform of a Poisson process with rate of 1. 
In the third case 3 = 0 , so 



m f (s) 



1 + s 



s + s 



again the Laplace transform of a Poisson process with a rate 
of 1. 

In the final case apply the Tauberian Theorem 



1+ (l+4S-2a6) s+(36+5 6 2 -2a6-5a6 2 +a 2 3 2 )s 2 



lim sm f (s) 
s+0 r 



2 3 222233 

+(26 +26 -3aB +a 8 +a 8 ) 



l+(l+46-3a6+a 2 3 2 ) s+ ( 36+56 2 -2a6-7a6 2 +3a 2 S 2 ) s 2 



2 3 222233 

+ (26 +26 -3a6 +a 6 +a 8 ) s 



1 

I' 



as required. 
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6 . The Spectrum of Counts 

For the statistical analysis of series of events the 
most useful quantity associated with a process is the 
(Bartlett) spectrum of counts. The spectrum of counts, 
g + (oj), is the Fourier transform of the covariance density of 
Nf(t). It is related to the Laplace transform of the inten- 



sity function, m^(s), by the relationship derived by Cox and 
Lewis [Ref. 29] 



We now derive this for the NEMA(l) process using III. C. 5. 4. 




— (1 + m, [iw] + rru [-iw] ) . 
it r r 



In that expression for m^(s), let 



a 



1 



1 + 48-2ctB 



(III.C.6 .1) 



38 + 58 2 - 2a8" 5aS 2 +a 2 6 2 



(III.C.6. 2) 



a 



3 



2 3 2 32223 

28 + 28 - 3a8 - 3a8 + a 8 + a 8 



(III .C. 6 .3) 



b 



1 



1 + 46 - 3a6 + a 2 8 2 



(III .C. 6 .4) 



b 



2 



38 + 58 2 - 2aB - 7a8 2 + 3a 2 8 2 



(III.C.6. 5) 



b 



3 



2 3 2 32223 

28 + 28 - 3a6 - 3aB + a 8 + a 8 . 



(III.C.6. 6) 



Then 
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, 2 3 

* 1+a, s+a~s +a.s 

m f (s) = ~ 2 - '3 4 

s+b^s +b 2 S +b^s 



Recall that X = 1 so 



2 3 

, 1+a, (iu) +a 9 (iu) +a.(iu) 

g + (a)) = i[l+ T 4 

iu+b (iu) +b (iu) J +b (iu) 4 



2 3 

1+a-^ (-iu ) +a 2 (-iu ) +a 3 (-iu) 

2 3 4 

-iu+b^ (-iu) +b 2 (-iu) +b 3 (-iu) 



] 



g. (u) = - 

3 + 7T 



[ioa+b 1 (iu) 2 +b 2 (iu ) 3 +b 3 (iu) 4 ] 

^ ^ *[-iu +b^ (-iu) 2+ t> 2 (-iu) 3+b 3 (-iu) 4 ] 



71 ' [iu+b, (iu) 2 +b„ (iu) 3 +b 0 (iu) 4 ] 



2 3 4 

(x [-iu+b^ (-iu) +b 9 (-iu) +b 3 (-iu) ] 



2 3 

[1+a^ (iu) +a 2 (iu) +a^(iu) ] 

+ x [-iu+b x (-iu) 2 +b 2 (-iu) 3 +b 3 (-iu) 4 ] 
[iu+b^ (iu) 2+ t> 2 (iu) 2+ t >3 (iu) 4 ] 

x [-iu+b 1 (-iu) 2 +b 2 (-iu) 3 +b 3 (-iu) 4 ] 



2 3 

[1+a^ (-iu) +a 2 (-iu) +a 3 (-iu) ] 

2 3 4 

x[iu+b 1 (iu) +b 2 (iu) +b 3 (iu) ] 

[iu+b 1 (iu) 2 +b 2 (iu) 3 +b 3 (iu) 4 ] 



x [-iu+b^ (-iu) 2 +b 2 (-iu) 3 +b 3 (-iu) 4 ] 



(III .C. 6 . 7) 
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Consider the first term with numerator and denominator the 



same. Let the denominator = D. 



D = [ iw+b-^ ( iu>) 2 +b2 ( ito) 2 +b 3 ( iw ) 4 ] [-ito+b^ (-ioo) ^+b^ ( - iw) ^+b^ ( - ico) ^ ] 



2 34 5 324 5 64 

to -ib 1 u) — b +ib^w +ib^w +b^oj -ib^b 0 w -b^b^to -b2to 



,. . . 5 , ,2 6 . , , 7 . , 5 , , 6..,, 7 , , 2 8 

+ib^b2<jJ +D2 a) “i^b^ -ib 3 u} -b-^b^oj +ib2b 3 w +b 3 to 



353 5 5 7 75 

l(-b-^u) +b 3 w +b^w -b^^w H-b^^w -b2b 3 u> +b2b 3 u> -b^w ) 



2 , 4 , , 2 4 , , 6 , 4,2 6 , , 6^, 2 8, 

+ (a) -b 2 u> +b^w -b^b^to - b2^ + ^ 2 t0 + ^ ) 3 a) ) 



a) 2 [1+ [b 2 -2b 2 ]a) 2 + [b 2 -2b^b 3 ]u) 4 +b 2 co 6 ) 



(III.C.6.8) 



Let 

X = l+(b 2 -2b 2 )a) 2 +(b 2 -2b 1 b 3 )a) 4 +b 2 a) 6 , 

where b n , b 2 > and b 3 are defined by III. C. 6. 4, III. C. 6. 5, 
III. C. 6. 6, respectively. Then 



D 



uj 2 X 



(III .C. 6 .9) 



Consider the numerator of the second term in 1 1 1. C. 6. 7 and 
call it N2 . Then 
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N2 = [1+a^ (ito) +a 3 (ito) 2 +a 3 (ito) 2 ] [-ioo+b^ (-iw) 2 + b 3 (-iu3) 2 +b, (-ioo) 



2.342 3 4 5 3 

-ico-b^u +0^2^ +b 3 o) +a^co -ia^b^co -a^b2W + ia^b 3 o) +ia2W 

+a2b^(jo^-ia2b2W^-a2b2U^-a2a)^+ia2b^a)^+a2b2a)^-ia2b2^^ 



N2 = i (-0)+ [a2~a^b^+b2 ] w 2 + (a^b2~a2b2+a2b^) w^-a^b^u^) (III. C. 6. 10) 

2 4 6 

+ ( a i~ b i ) a) + (a 2 b 1 -a 1 b 2 -a 3 +b 3 ) oj + (a 3 b 2 -a 2 b 3 ) to , 



where a^^ , a 2 , a 3 , b 1 , b 2 , and b 3 are defined in III. C. 6.1 
through III. C. 6. 6 respectively. Consider the numerator of 
the third term in III. C. 6. 7 and call it N3. Then 



N3 = [1+a-j^ (-iu) +a 3 (~ico) 2 +a 3 (-iu) 2 ] [ico+b^ (iw) 2 +b 3 (ico) 2 +b 3 (iu) 



• u 2 -u 3^, 4^ 2 . 3 , 4 . , 5 . 3 

lco-b^o) -ib2to +b 3 0J +a^co +ia^o^03 -a^b2^ -la^b^ -ia20) 



, . 4 , * , 5 ,6 4 . . 5. . 6.. . 7 

+a o b,0) +ia~b~u3 -a_b_,co -a^co -ia-b-,03 +a-,b-o) +ia_,b o (j0 
21 22 23 j 31 32 33 



3 5 7 

N3 = i (co- [ a 3 — a^b ^+b 3 ] w - [a^b 3 ~a2b2+a 3 b]_ ] co ta^b^ ) (III. C. 6. 11) 



2 4 6 

+ (a^-b^)co + (a 2 b^-a^b 2 ~ a 2 + b 3 ) co + ) 60 , 



where a^ , a 3 , a 3 , b-^, b 3 , and b 3 are defined in III. C. 6.1 
through III. C. 6. 6, respectively. Note from III. C. 6. 7 that 
all terms in the sum have the same denominator. Use III. C. 6. 10 
and III.C.6.11 to determine the numerator of the second and 
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third terms and call it N. Then 



3 5 7 3 

N = i [-oj+ ( a 2 -a ]_k^+b 2 ) w + (a^b 2 -a 2 b 2 +a 2 b^ ) co -a^b^co + co- (a 2 -a^b^+b 2 ) co 

5 7 2 4 

- (a^b^-a^^a^b-^) co ta^b^co ] + (a^-b^)oo + (a 2 b 1 -a^b 2 -a 2 +b 2 ) co 

+ (a 2 b 2 ~a 2 b 2 ) oo^+ (a^-b-^) to 2 + (a^-^-a^b^a^+b^ ) co 4 + (a 2 b 2 ~a 2 b 2 ) to® 

= 2 [ (a^-b 1 ) oo 2 + (a 2 b^-a^b 2 -a 2 +b 2 ) co 4 + (a^b^a^^) co 6 ] . 



Let 

2 6 

y = ( a ]_”bj_) + ( a 2 b^-a^b 2 ~a 2 +b 2 ) co + (a^b^a^^) co ] . (III. C. 6. 12) 

Then 



N = 2co y 



(III. C. 6 .13) 



Using III. C. 6. 7, III. C. 6. 8, III. C. 6. 13 



g, (u) = 



,2 0 2 

1 /Co x 2co 

7T 1 2 2 

CO X CO X 






= I(^), 



IT X 



(III .C.6 .14) 



where x and y are defined in III. C.6. 9 and III. C.6. 12, 
respectively. 
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Figures III. C. 6.1 through III. C. 6. 3 show the results 
of the calculation of the (Bartlett) spectrum of counts. In 
presenting the results the constant ^ in III. C. 6. 14 was ig- 
nored. Figure III. C. 6.1 shows the spectrum of counts for the 
same a and 6 values that were used for the sample paths and 
scatter plots of Figures III. C. 3.1 through IIIC.3.6. This 
figure also shows the variation in the spectrum of counts as 
the P > X^) varies from its lowest to highest values. 
Figure III. C. 3. 2 holds the P(X . > x ) constant and varies 
the correlation. Since the spectrum of counts for a Poisson 
process is a constant one when A equals 1 and the constant 
— is ignored, the correlation can be viewed as a measure of 

TT 

the process' departure from a Poisson process. This 

divergence as a function of the correlation shows clearly in 

this figure. Figure III. C. 6. 3 holds the correlation constant 

and varies the p ( x n+ j_ > x n ) • The slight variation in the 

spectra shows that while the spectrum of counts does vary 

with the P(X ,. > X ) , the correlation plays a more dominant 
n+1 n 

role . 

The analysis from the Laplace transform of sums in 
III.C.4, through the Laplace transform of the intensity func- 
tion in III.C.5, to the spectrum of counts in this section 

1 

can be performed using the correlated {E^} , {E n } sequences 
of III.C.2 and thus for negative correlations. Details are 
not given. 
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FIGURE III. 
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FIGURE III. 
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FIGURE III. 



7 . Joint Laplace-Stieltjes Transform of X ^ and X ^^ 



Because the NEMA(l) process is, by construction, 
only one-dependent, all of the second-order properties of 
{X n > are contained in adjacent pairs {X n ,X n+ ^}. previous 

sections quantifiers of the distribution of {X n ,X n+ ^} such 
as p. and P(X , . > X ) have been derived. Here we give the 
Laplace-Stieltjes transform of the joint distribution. One 
could, for example, study the effect of the two parameters 
from this result by deriving directional moments. 

The joint Laplace-Stieltjes transform of X n and X n+ ^ 
can be calculated by considering each of the sixteen possible 
combinations of parameter values for X n and x n+ j_/ as was done 

1 _ g 

in III.C.3. Let - r - r = 6, (1-a) 0 = y, cj> x x (s.,s 2 ) = 

-s.x -s-x . . n ' n+1 

E(e n n X ) , and <J> E (s) = Then 
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<J> X x ( s 1 ' s 2 ) = aa66E(e 

n' n+1 



- s l lE n +SE n-l ! - s 2 lE n + l +6E n J 



-s, [E +6E , ] -s„E , 

+ (l-a)a66E(e 1 n n 1 2 n+1 ) 

-s, E -s 9 [E , +8E ] 

+ a (1-a) 55E (e ) 

-S-.E -S-E , 

+ (1-a) (1-a) 66E(e 1 n z n x ) 

-s, [yE +8E , ] -s_ [E , +8E ] 

+ aa6 (l-o ) E (e ) 

-s, [yE +BE , ]-s_E ,, 

+ ( 1-a ) a6 (1-6 ) E (e 1 n n_1 2 n+1 ) 

-s,yE -s 9 [E -,+BE ] 

+ a (1-a) 5 (1-6 ) E (e 1 n 2 n+1 n ) 

-s,yE -s~E , 

+ (1-a) (1-a) 6 (l-5)E(e 1 n z n ^ x ) 

-s, [E n +6E , ]-s 9 [yE , +3E ] 

+ aa (1-6 )oE(e ) 

-s, [E +6E ,1-s-yE , 

+ (1-a) a (1-6) 5E(e 1 n n_1 2 n - 1 ) 

-s, E -s_ [yE , +BE ] 

+ a (1-a) (1-6 ) 5E (e 1 n 2 1 n ) 

-s,E -s yE , 

+ (1-a) (1-a) (1-6) 6E(e n 2 ) 



+ aa ( 1-6 ) (1-6 ) E (e 



-s 1 [yE n+ 8E n _ 1 ]-s 2 [yE n+1+ 3E n ]^ 



+ (1-a) a (1-5) ( 1-6 ) E (e 



-s 1 [yE n +eE n . 1 ]-s 2 yE n+li 



+ a ( 1-a ) (1-6) (1-6 ) E (e 
+ (1-a) (1-a) (1-6) (1-5 ) E (e 



- S l YE n- S 2 [YE n-l +6E n ] 



- s l yE rT s 2 Y W 
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aa6 6<}> E (s^Bs^ e E (Bs^) e £ ( s 2 ) 



+ (l-a)a66 4) E (s 1 ) (J) e (6s 1 ) <p E (s 2 ) 

+ a (1-a) 66 4> E (s 1 +Ss 2 ) <{> E (s 2 ) 

+ (1-a) (1-a) 66<J> E (s^) <{> E (s 2 ) 

+ aa6 (1-6 ) 4 > e (ys 1 +6s 2 ) <J> E ( 3 s 1 ) <|> E (s 2 ) 

+ (1-a) a6 (1-6 ) <() E (ys 1 ) <j> E (Ss 1 ) <j> E (s 2 ) 

+ a (1-a) 6 (1-6 ) <^ E (ys 1 +6s 2 ) ()) E (s 2 ) 

+ (1-a) (1-a) 6 (1-6) 4 > e (ys 1 ) 4> e (s 2 ) 

+ aa (1-6 ) 6 4> E (s^Bs^ <p E ( 3 s ^ ) <j> E (ys 2 ) 

+ (1-a) a (1-6 ) 6 <j> E (s 1 ) <j) E (Bs 1 ) <J> E (ys 2 ) 

+ a (1-a) (1-6) 6 <|> e (s^+6s 2 ) <^ e (ys 2 ) 

+ (1-a) (1-a) (1-6)6 <Ji e (s 1 )<|) e (ys 2 ) 

+ aa ( 1-6 ) (1-6) <|> e (ys 1 +6s 2 ) 4> e (Bs 1 ) 4» e (ys 2 ) 
+ ( 1-a) a ( 1-6 ) (1-6) <(> e (ys 1 ) 4> e (Bs 1 ) <() e (ys 2 ) 
+ a (1-a) (1-6) (1-6) <j> E (ys^Bs^ <t> E ( Y s 2 ) 

+ (1-a) (1-a) (1-6) (l-6)<|» E (YS 1 )(f» E (YS 2 ) 
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[6 (f> E ( s 2 ) h- ( 1-5 ) $ E (ys 2 ) ] 



x [aa6<|) E (s 1 +gs 2 ) (J E (gs 1 ) + (1-a) a6 $ E (s 1 ) 4> E (gs.^) 



+ a (1-a) 5(J) E (s 1 +gs 2 ) + (1-a) (1-a) 6<f> E (s^ 



+ aa (1-6 ) 4> E (ys 1 +gs 2 ) <j> E (gs 1 ) 



+ (l-a)a (1-5) <J> E (ys 1 ) <j> E (gs 1 ) + a (1-a) (1-6 ) <J> E (ys^+gs^ 



+ (1-a) (l-a)(l-6)4) E (ys 1 ) ] 




[64> e (s 2 ) + (1-6) <t> E (y s 2 ) ] 



(III. C. 7.1) 



x [a<j) E (gs) + (1— a) ] x [ao(j) E (s 1 +gs 2 ) + (1-a) S<f> E (s^ 

+ a(l-5)(J> E (ys 1 +gs 2 ) + (l-a) (1-5 ) <j> E (ys 1 ) ] 

For the special cases of the EMA(l) process, III. C. 7.1 reduces 
to the results given in Lawrence and Lewis [Ref. 5]. 

D. THE MOVING MINIMUM MODEL 
1. Introduction 



one-dependent sequences of random variables with marginal 
Exponential distribution is the so-called minimum model. With 
this model the {X n } sequence is generated by taking the moving 
minimum value of two Exponential random variables . The proposed 
generation scheme is 



Another possible scheme that can be used to generate 




MIN (E ,bE . 
n n— 1 



(III .D. 1.1) 
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where (X^, n = 1,2,...} is a sequence of random variables 
with marginal Exponential distribution, (E n , n = 0,1,...} is 
an iid sequence of Exponential random variables with unit 

mean, and b >_ 0 . This will produce an {X^} sequence with a 

Id 1 Id 

rate of — and an expected value of This expected 

value produces one difficulty since it is a function of the 

parameter, b. This complicates comparisons between results 

with different parameter values and decreases the value of 

scatter plots and sample paths. However, this difficulty 

can be easily removed by multiplying the {x^} by — jj- . The 

generation scheme then becomes 

X n = X' = MIN( [^]E , [b+l]E ,), (III. D. 1.2) 

non bn n- 1 

with (E } and b defined as before. The (X } has a rate of 
n n 

one and, hence, an expected value of one. This facilitates 
comparisons for different parameter values with the NEMA(l) 
discussed in III.C which produces random variables with unit 
means . 

The investigation of the moving minimum model is moti- 
vated by the previous result in III.C. 2 that linear additive 
models have a constrained range of serial correlation. The 
hope is that the non-linearity of the moving minimum model 
will obviate this constraint. The minimum scheme has been 
used by Tavares [Refs. 22 and 30] to generate first-order 
autoregressive exponential processes and by Marshall and 
Olkin [Ref. 31] to generate correlated bivariate Exponential 
variables . 
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2. Correlation Structure 



The first-order serial correlation can be computed 
by the following approach 



E(X n+l X n ) 



E ( [MIN ( E n+1 , [b+l]E n ) ] 

[MIN ( t^]E n , [b+l]E n _ 1 ) ] ) . 



The terms inside the expected value can be made independent 
by conditioning on the value of E n . The E(X n+ .^X n ) is then 
found by multiplying the conditional result by the density 
of E , and integrating. Implementing this approach we have 



E(X n+l X n ) = / E([MIN([^]E n+1/ [b+l]y)] 



x [MIN ( [^!p]y , [b + l]E n _ 1 ] | E n _ ) e y dy 



E(X ni .X) = / (E [MIN ( [— rj— ] E . , [b + 1] y ) ] ) 



n+1 n 



b J n+1' 



(III .D. 2. 1) 



r b+l 



x(E[MIN( [~i]y, [b + l]E n _ 1 )])e Y dy 



The expected value of the minima can be calculated as follows 



E(MIN[F±i)E ,(b+l)y]) = / x (x^r) e “ ’ ^dx 



,, , . , bx 

(b+1) y _ „ 

b+1' 



0 



bx 



+ / 



(b+1) y 



(b+1)y( bTl )e b+ldx 
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bx 



E (MIN [ (~£~) E n+1 / (b+ 1 ) y ] ) = -xe b+1 



(b+1) y 



0 



(b+1) y . 

+ I e dx + (b+1) ye ^ 

0 



-by A ^W^ b+1)y (^ )e 'b+ X l dx 



= - (b+1) ye -b -^ + (S+l)/ 



+ (b+1) ye b ^ 



E (MIN [ (£±I) E n+1 , (b+ 1 ) y ] ) 



— (1 - e -b ^) 



(III .D.2. 2) 



Similarly, 



E(MlN[(Sii)y, (btllE^jl) 



(b+1) y/b x 
/ x( bTT )e dx 



+ / ib+ili, 1 )e b*l dx 

(b+1) y/b b b+1 



= -xe 



b+1 



x 



(b+1) y/b (b+1) y/b -5+3- 

+ / e dx 

0 0 



(b+1) y -y/b 
b 



= -ib+llx g-y/b + (b+1) y 



X 



(b+1) y/b x - 5TT 



( b+T } e dx 



(b+1) y -y/b 
b 
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(III. D. 2. 3) 



E(MIN[ (b+l)E n _ 1 ]) = (b+1) ( 1 - e _y/b ) . 



Using III. D. 2. 2 and III. D. 2. 3 in III. D. 2.1 produces 



E(X n+l X n> 



= / <!£!> (b+1) (l-e-y /b )e -ydy 



(b+1) b+1 



— (b+1) + 



(b+1) 



b ( 1+b+g-) 



(b+1) 

b 2 +b+l 



Therefore, since E(X ) = 1 

n 



COV(X n , 1( X ) = 1 

n+1 n b 2 +b+l 



b 2 +b+l 



and 



CORR(X . ,X ) = (III. D. 2. 4) 

n+i n b Z +b+l 

Thus the model allows a range of correlations from 
[0,-j] . The minimum value is achieved when b is zero or in 
the limit as b tends to infinity. The maximum value is achieved 
when b is equal to one. An interesting aspect of the correla- 
tion structure of the moving minimum model is that reciprocal 
values of b produce equal correlations. This is a similar 
kind of "invertibility " found for the other moving average 
models discussed in III.C. The range of b could be restricted 
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to the interval [0,1] without reducing the possible range of 
correlations. However, doing so, as with the NEMA(l) model, 
would ignore the fact that in the non-normal case character- 
istics other than correlation may be different in the case of 
b and Also, as is the case for all the first-order moving 

average processes addressed in this paper, the correlation 
for lags greater than one are zero. So the study of correla- 
tion structure is limited to the study of the serial correla- 
tion with lag one. 

3 . Negative Correlation 

The range of possible correlations can be extended in 
a fashion similar to the NEMA(l) model (see III.C.2) by the 
use of correlated or antithetic variables. Using this approach 
the generation formula becomes for antithetic variables 

X n = MIN( [^jji]E n , [b+l]E^_ 1 ) , (III. D. 3.1) 

where all variables are defined as in III. D. 1.1 and {E^, 

n = 0,1,...} is generated from the {E } sequence using the 

-E n 

relationship E^ = - In (1-e n ) . Note that this implies that 
{E^} is also iid Exponential with unit mean. Again 

E(X n x n+ i ) = E[(MIN[(^i)E n , (b+DE^J) 

x (MIN[ (^~)E n+1 , (b+l)E^] ) ] 

and conditioning on the value of E^, multiplying by the density 
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of E^, and integrating produces 



E(X n X n+l ) = / BdMINd^ly.Ib+llE'^)] 



rb + 1 



[MIN ( [i^±] E n+1 , [b+1] E^) ] | E n = y ) e y dy 



, K + 1 

E(X n X n + l ) - / q (E [MIN ( [=£±]y, [b+l]E^_ 1 ) ] ) 



(III. D. 3 .2) 



rb+1 



(E [MIN ( [^±] E n+1 f - [b+1] In [l-e _y ] ) ] ) e~ y dy 



The first expected value is identical to III. D. 2. 3. Thus 



E(MIN[ (^)y, (b+l)E^_ 1 ] ) = (b+1) (l-e" Y/b ) 



(III. D. 3. 3) 



The second can be calculated as before 



E (MIN [ (~£~) E n+1 /“ (b+1) In ( 1-e Y ) ] ) 



- (b+1) In ( l-e -y ) fa 
= / Q x( b¥I )e dx 



bx 

OO “ — =- 

+ j -(b+l)ln(l-e y ) e + ax 

- (b+1) In (l-e _Y ) 



= -xe 



bx 

‘b+1 



- (b+1) In ( l-e -y ) 



-(fa+Uind-e-l'l 

( r ) / (tttt) e ax 



0 



b+1 



- (b + l)ln(l-e-y)e bln(1 ‘ e ’ y) 
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E (MIN [ (— j“) E n+1 / - (b+1 ) In ( 1-e y ) ] ) 



(b + l)ln(l-e- y )e bln(1 - e " Y) + (^ ( l-e bln [1 - e ~ Y] ) 



- (b + l)ln(l-e- y )e bln(1 - e " Y) 



E(MIN[(H±i)E n+1 ,- (b+1) In (l-e _Y )]) 

= (^p) (1 - [l-e _y ] b ) (III. D. 3. 4) 

Substituting III. D. 3. 3 and III. D. 3. 4 into III. D. 3. 2 yields 



E(X n X n + l> 



= /"<fe+i> U-e' y/b ) (6ji) (l-[l-e-i'j b )e -l'dy 



_ b+1 

,b+l^ 2 ,b+l N r°° ,b+l x b Y j 

( ~b~ ) I e dy_ { ^b~ ) J { ~b~ )e dy 

u 0 0 



-b+1. 2 ( -y -y. b , 

- (-T— ) } e J (1-e T ) dy 

D 0 



,b+l,2," -< 1+ E>y -y, b , 

+ (~u~ ) J e (1_e > dy 



The first two integrals are trivial. In the third the change 
of variable z = (l-e~ y ) , dz = e y dy makes that integral straight- 
forward. In the last integral, the change of variable u = e y , 

— — = dy makes that integral recognizable as the integral of 
a Beta random variable. Using these changes of variables and 
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making the appropriate changes to the limits of integration 
produces 



E (X X « ) 
n n+1 



Then 



(b+1) 2 _ b+1 



- <^) 2 1 



b , 

z dz 



0 -(l+£) . , 

/ u b (l-u) b (-^) 

b u 



,b+l, ,b+l, b+1 , c 1 1/b,, 

(- 5 T— ) - (-C— ) y + J u ( 1 — u) du 

D D b z 0 

r (i+i) r (i+b) 

r (2+b+g.) 

g-r (g)br (b) 

( 1 +b+i) (b+i) r (b+i) 

b 2 r (i) r (b) 

(b 2 +b+l) (b 2 +l) r (b+i) 

D 



COV(X ^ X . n ) 
n' n+1 



b 2 r (g-) r (b) 



(b 2 +b+l) (b 2 +l) r(b+i) 



- 1 



and 



C0RR(X n' X n+l ) 



b 2 r (g-) r (b) 

(b 2 +b+l) (b 2 +l) r (b+g-) 



(III .0.3 .5) 



Like the expression for positive correlation, this 
expression is also symmetric with respect to reciprocal values 
of the parameter. It attains a minimum value of minus one- 
third when the parameter value is one. Graphs of the 
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correlation as a function of the parameter b for both posi- 
tive and negative correlations are provided in Figures III. D. 3.1 
and III. D. 3. 2, respectively. Unlike the NEMA(l) model which 
requires an additional parameter (see equation III.B.6) to 
achieve a full range of negative correlations, the moving 
minimum model can achieve its full range with the single 
parameter b and an antithetic sequence {E^}. 

4 . Joint Density of X ^ and X n ^ 



possible using a conditioning argument to determine 

P ( X_ < x|E = z) and P(X .. < y|E = z) . These values along with 
n — 1 n n+1 — 1 n 

the probability that E^ takes on a given range of values are 
sufficient to determine the joint distribution function of X 
and x n+ ]_* T he form of the distribution will vary depending on 
whether one is above or below the line X . = bX . The joint 



density, where it exists, is determined by differentiating 
the distribution function. 



Calculation of the joint density of X and X , . is 

n n+1 



From III. D. 1.2 we have 



X. 



n 



MIN ( [— g— ] , [b+l]E n _ 1 



Then 




1 




P (X < x| E = z) 
n — 1 n 



(III .D.4 .1) 



V 1-e b+1 if (^Ji) z > x 

b 
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.00 1.25 2.50 3.75 5.00 6.25 7.50 8.75 

PRRAMETER VALUE 



The first result in the above is obvious. To justify the 
second, consider 



P(X n <x| E n=z>^-) 



P ( [b+1 ] E . < x) 
n-1 — 



- P(E n-l -bfr’ 

X 



= 1 - e 



b+1 



Since X n+1 = MIN( [^]E n+1 , [b+l]E n ) , then 



P(X . < y E = z) = 
n+1 — 1 1 n 



. I 



if z < 



. b+1 . - 

1-e if z > 



b+1' 



b+l‘ 



(III. D. 4 .2) 



To consider the joint distribution, note that when 

(i.e. when you are above the line bx = y) , the range 
of possible z values can be broken up into three regions. See 
Figure III. D. 4.1. Then, since E^ is Exponential with unit 
mean. 



P(z 6 REGION 1) 



P(z e REGION 1) 



P (z < 



bx 



-b+1' 



1 



bx 

b+1 



(III .D .4 . 3) 



P(z e REGION 2) 



P ( 



bx 

b+1 




f 
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FIGURE III. 



P(z € REGION 2) 



P(Z £ REGION 3) 



P(z e REGION 3) 





) 



Now by definition 



(III. D. 4 .4) 



(III .D .4 .5) 



p ( x n 1 X ' X n +l £ Y I 2 e REGION i) = P(X n £x|z £ REGION i) 

X P (X n+l- y l Z e REGI0N i) (III. C. 4. 6) 

because when conditioned on the value of E , these probabili- 
ties are independent. Using the above equation, III. D. 4.1, 

III. D. 4. 2, and the definition of the regions in Figure III. D. 4.1, 



P(X < x, X .. < y I z £ REGION 1) = 1 (III. D. 4. 7) 

n — n+i — 

x 

P(X n < x,X n+1 < y| z £ REGION 2) = 1 - e b+1 (III. D. 4. 8) 

£_ 

P(X n < x,X n+1 < y|z £ REGION 3) = (1-e b+1 ) (1-e D+1 ) (III. D. 4. 9) 

Using the results of III. D. 4. 3 through III. D. 4. 5 and III. D. 4. 7 
through III. D. 4. 9 we can compute the joint distribution of X n 
and when y > bx by using the relation 
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(III. D. 4 .10) 



P ( X n ix, X n +1 <y> - 


3 

l P (X < x, X , < y | z e REGION i) (III.D.4 .10) 

1=1 n n 1 

x P (z e REGION i) 

bx x bx 

1 (l-e b+1 ) + (l-e b+1 )(e b+1 . e -y/bH-l, 

+ (l-e b+1 ) (l-e b+l )e -y/(b+l) 

bx bx y (x+y) y 

, b+1 , b+1 -x b+1 , b+1 ^ 5+1 

l-e +e -e -e +e +e 

-rrf -y -< b }i>(* + (b+i]y> 

- e -e -*+e 


P(X ni K 'Vlif - 


X 

l-e“ x -e" Y +e b+1 Y (III.D.4. 11) 



Similarly, when y < bx (i.e. when you are below the 
line bx = y) , the range of possible z values can be separated 



into three regions . 


See Figure III.D.4. 2. Then 


P ( z s REGION 1) 


= p(2 ibTr>' 


P ( z £ REGION 2) 


- l-e b+1 , (III.D.4. 12) 


P(z £ REGION 2) 


- E>< b$r < 2 -b+i 1 ' 


P(z £ REGION 2) 


y _ bx 

b + 1 b+1 / TTT t~i 4 1-3 1 

= e - e . (III.D.4. 13) 


P(z £ REGION 3) 


= P < 2 > b+ X l } 
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P(z € REGION 3) 




(III .D. 4 .14) 



Using III. D. 4.1, III. D. 4. 2, III. D. 4. 6, and the definitions of 
the regions in Figure III. D. 4. 2, the following results hold 
for the given region. 



P (X <_ *,X n+1 £y | z s REGION 1) = 

P (X„ < x,X ^ < y I z e REGION 2) = 
n — n+l — 1 

P(X <x,X . <y|z e REGION 3) = 

n — n+l — ■* 1 



X 



x 



.. b+1. 

(1-e ) 



1 




, , b+1. 

(1-e ) 



(III .D.4 .15) 
(III. D. 4 .16) 
(III .D.4 .17) 



Combining III. D.4. 10 with III. D.4. 12 through III. D.4. 17 yields 
for bx < y 



P (X <X/X .. <y) 
n — ' n+l — 1 



-x-^L 
. -y -x b+1 

i-e - 1 -e +e 



(III. D.4 .18) 



x 



Let F (x,y) = 1-e -e y +e , the distribution 



function of X and X when bx < y; and let f (x,y) be the 
n n+l 

joint density of X n and X n+ ^ when bx < y . Then 



f (x,y) = 



kl7 Fl(x 'yi ■ 

X 

3 . -y b+1 ' , 

K le e 1 



X 



3 3 ri -x -y, b+1 1 . 

[1-e -e “ +e ] 

3x 3y 
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f 1 (x,y) 



bx < y; 



x > 0 . 



(III. D. 4 .19) 



, 1 , b+1 y 

= { b+T )e 



For bx ,< y, the distribution function, F (x,y) , is 

1 . e -y. e - +e - x -Btr .,,2 



If f (x,y) is the joint density of X 



n 



and X n+ -^ when bx > y, then 



f (x,y) = 



3 3 „2 , . 3 3 ri - y “X , 

r- 7 - F (x,y) = tt— -k— ll - e J -e +e 

3x 3y '■* 3x 3y 



-x-^L 



b+l ■ 



f (x,y) = 



- X _^L 

— [e - ^- (-^— ) e b+1 ] 

3x ie b+1 ' e J 

-x-M- 



( bTI )e 



b+1 



y < bx ; y > 0 . 



(III.D. 



Note that there is a positive probability that the 
point ( x n ' x n+ ]_) lies on the line bx = y . This probability 
can be computed as follows. We have 



X n = MINI hrlE n , [b+1 ]£„_!> 

X n+1 “ MIN( l3 i l E n + l' [b+1IE n ) 



The point (X ,X , . ) lies on the line bx = y when X = ( . — ) E 
r n n+1 n b 

and X , . = (b+1 ) E . Now 
n+1 n 



P(X = [— ^]E ;X .. = [b+1 ] E ) 
n on n+1 n 



- p ”X 1E ni IbtllE n-r ‘ b+11E n i >T 1E n + l> 



4.20) 
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The events in the right hand side can be made independent by 
conditioning on the value of E n « Then 



,b+l 



P(X n= ]E n ;X n + l= t b+11 V 



= ( [ T ] y - tb+1 1 E n-1 ' [b+1 1 y - [ T ] E n+1 1 E n = Y * e Ydy 



• b+1 ■ 



= / 0 P(E n-l > b )P(E n+l >b y )e “ Yd y 



^"1 / (b+l+i) e 



l - (b+l+g-) y 



dy 



b+l+^ 0 



P(X n = 'Tr>V X n + l= ' b+11E n> 



b 2 +b+l 



(III.D.4 



Because there is a positive probability of lying on the line 
bx = y, the moving minimum model can be said to have a line 
degeneracy. An important implication of the positive proba- 
bility of ( x n ' x n+ ]_) lying on the line bx = y is that the 
moving minimum model will produce runs of values of constant 
ratio b. The values of (X n ) in these runs will be decreasing, 
equal, or increasing for b less than, equal to, or greater 
than one, respectively. The length of the runs will be geo- 
metrically distributed with parameter „ b for the positive 

b +b+l 

correlation case. It was this kind of degeneracy in the 
Exponential autoregressive model, EAR ( 1) , that proved to be 
one of the model's major weaknesses. The degeneracy also 



. 21 ) 
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occurs in the Tavares autoregressive model and the bivariate 
Exponential pairs derived by Marshall and Olkin. 

The probability of lying above or below the line bx = y 
can be easily found by integrating the appropriate joint den- 
sity over the area desired. Thus, for bx < y, 



P (lying above bx = y) = 



I / f 1 (x,y)dydx 
0 bx 



x 



( ( t 1 \ b+1 

/ / ^ e 



0 bx 



dydx 



r , 1 > b+1 -bx, 

/ <S7T )e e dx 



P (lying above bx = y) = 



b +b+l 



(III. D. 4. 22) 



Similarly for bx > v 



00 bx 



P (lying below bx = y) = / / r (x,y)dydx 

0 0 



CO bx 



= ; o ( ^ )e 



by 

b s _ -x 5+T 



dydx 



b 2 x 

= / e _x (1-e b+1 ) dx 
0 



= 1 - 



b+1 



b +b+l 
. 2 



P (lying below bx = y) = 



b +b+l 



(III .D.4 .23) 
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5. Conditional Expectation and P(X > X ) 

n+1 n 

Besides the correlation coefficient, there are two 
other characterizations of the joint distribution of X and 
X n+ -^ which we have considered. They are the conditional 
expectations and the P (X ^ > X n ) . Both of these quantities 
can be derived by considering the four possible sets of values 
for X n and X n+ ^ , computing the probability of each set occurring, 
and weighting the conditional expectation or probability by 
its probability of occurrence. 



First, the probability of occurrence for each set of 



values must be calculated. Consider the case where X 



n 




and X 






By conditioning on the value of E , the events on the RHS 
become independent. The calculation then proceeds in a 



straightforward way . 




/ o P([ ^ ] Y£ [b+l]E n _ 1 , [^]E n+1 < [b+l]y|E n = y)e _Y dy 




0 



n-1 b 



> J: ) p ( E n+1 £ by ) e“ Y dy 
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P(X n = [^]En.X n+1 = I^)E n+1 ) 



“ -(r^Dy 
/ e dy 

0 



» -(b+l+±)y 
/ e b dy 

0 



b+1 b 2 +b+l 



!(I n =l T 1! .'Vl- [ T 1 W 



(b+1) (b^+b+1 ) 



(III. D. 5.1) 



The second case when X = (^-^) E and X , . = (b+l)E has 

n d n n+l n 

already been computed. The result is III. D. 4. 21 and is re- 
peated here as 



r b+l 



P(X n= : IE n'Vr !b+11E n> = 72 



(III. D. 5. 2) 



b +b+l 



L I I 

The third case is when X = (b+l)E . and X ,. = ( °, --) E . 

n n-1 n+l b n+l 

Here we proceed as in the first case. 



P(X n = (b + l)E n . 1 ,X n = I^]E n+1 ) 



P ( [b+1] E^ < [^±1] E n , [fe±i] E n+1 < [b + l] £„) 



/ Pdb+llE^ < [^]y, (^]E n+ il [b+l]y|E n = y)e y dy 



= / P(E n _ 1 < y )P(E n+1 < by)e" Y dy 
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b+1 



P(X n = tb + l]E n . 1 ,X n = [VlVl> 



oo -(l+_)y 

/ e y dy - J e 
0 0 



-Y ^ -f - ' b ' J dy -/ e _(b+1)y dy + J 



» -(b+l+±)y 



0 



dy 



1 _ _ 4 ^ + 1 



b+1 b + l b 2 +b+1 



b+1 



P(X n =[b+llE n-l' X n + f'T 1 Vl ) = 



b 2 +b+l 



(III.D.5.3) 



The final case is when X = (b+l)E , and X ,, = (b+l)E . As 

n n-1 n+1 n 

before 



P(x n = [b + l]E n .l' x n + l= [b+1]E n) 



/ P ( (b+1 ] E n _ 1 < l^±^]y, [b+1 ] y <_ l ii 5 i ) E n+ i l E n = y)e Y <Sy 



■ b+1 



/ 0 P < E n-lig )P,E n+l >b i’ )e ' ydY 



00 00 ~ (b+ 1 + 7 -) y 

/ e y dy - / e dy 



b+1 b 2 +b+i 



p (X n = tb+l]Vl' x n+l“ tb+l]E n ) 



(b+1) (b +b+l) 



(III. D. 5. 4) 



The conditional expectations can now be written by 
inspection. 
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Case 



X = [^i]E ,X , 
n bn n+1 



rktilE 

L b J n+1 



r b+l 



X n = ]E n' X n + l = ^ b+1]E n 



X n * lb+11E n-l' X n + l= [ ^ IE n + l 



X n “ [b+1JE n-l' X „ + l= ‘ b+1 l E n-l 



E(X n+1 lX n ^li 

b+1 

b 

by 

b+1 

b 

b+1 



Weighting these conditional expectations by the probabilities 
in III. D. 5.1 through III. D. 5. 4 yields the final result. 



E(X n+ll X n = y) = I E < X n+i I x r, = y;case i)P(case i) 



i=l 



n+1 1 n 



= [^3 [ ~2 ] + (by) (“^ ) 

D (b+1) (b z +b+l) b +b+l 



+ (^) ( ) + (b+1) ( K: ) 



b b 2 +b+l 



[b+1] [b +b+l ] 



E(X n+ll X n = y) = 



., b 2 y+l 

b 2 +b+l 



(III. D. 5. 5) 



It is quite surprising that the regression of X n+ -^ on X n is 
linear in y, considering the non-linearity of the process which 
generates the (X n ) . 

The conditional expectation of X n given X n+ -^ can be 
derived with equal dispatch. 
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Case 



£(X n iX n+1 ^d. 



x n- 



[ ] e 

L b J n+1 



x n = 'TT 1E n' x n + l 



= [b+1 ] E 



n 



x n = ' b+1 > E n-l' x n + l = 'T lE n rt 



x n = [b+11E n-l' x n + l= [b+1JE n 



b+1 

b 

Z 

b 

b+1 

b+1 



Using III. D. 5.1 through III. D. 5. 4 as before 



E(X n |X 



n+1 



y) 



4 

l E(X n |X n+1 = y,case i)P(case i) 
i=l 



- (■ 



[b+1] [b +b+l ] 



) - <£> (-j 6 —) 

D o +b+l 



+ (b+1) (- 



-) + (b+1) (- 



b +b+l 



[b+1] [b +b+l ] 



-) 



E(X nl X n+l 



= y) = 



2 

b +y 
b 2 +b+l 



+ 1 



(III.D.5.6) 



The probability that X n+ ^ is greater than X n can also 
be easily computed if one is careful to differentiate between 
the case where b < 1 and b > 1 . 



Case 



P (X 



n+1- 



■V 



x = [— ] E ,X . = [^^]E . 

n b J n n+1 b J n+1 



1 

2 
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X n= IE n' X n + l= lb+11E n 
X n = |b+11E n-l' X n+l= ^ !E n+l 

x n= (b+1 ] E n _i , X n+ ^ = [b+l] E n 

Thus when b £ 1 we have, using III. D. 5.1 



j 0 if b < 1, 
(l if b > 1 

1 

b+1 

1 

2 

through III. D. 5. 4 



P(X > X ) 
n +1 n 



i=l 



P (X . > X |case i)P(case i) 
n+l n 



<§>< 



[b+1 ] [b +b+l ] 



-) + <ctt> ( -v b ) 
b+1 b 2 +b+l 



+ (j) ( 



[b+1] [b +b+l] 



-) 



P{X n+l >X n } 



1 

2 



(b+1) (b +b+l) 



b < 1 



(III .D. 5. 7) 



A similar computation with b > 1 again using III. D. 5.1 
through III. D. 5. 4 yields 



P (X . , > X ) = i + , b > 1. (III. D. 5. 8) 

n+1 n 2 (b+1) (b 2 +b+l) 

Thus a graph of p ( x n+ ]_ > x p ) will have a discontinuity at b = 1 
when case 2 switches from a probability of zero to a probability 
of one. This graph is presented as Figure III. D. 5.1. The 
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FIGURE III. 




minimum value of one-third occurs at b = 1. The maximum value 
of two-thirds occurs at b = 1 + . The moving minimum model, 
therefore, has a greater range of values for the P (X n+ ^ > X n ) 
than does the NEMA(l) model. However, the greater range for 
the ? (X n+ ^ > X n ) and greater range of correlations must be 
balanced against the degeneracy of the model. 

As was noted with the NEMA(l) model, the correlation 
in non-normal models does not define the joint properties of 
X n and X n+1 « Although the cases of b and ^ are indistinguish- 
able from the viewpoint of correlation (see III. D. 2. 4 and 
III. D. 3. 5), these cases will have significantly different 
sample paths as indicated by III. D. 5. 7, III. D. 5. 8, and the 
discussion of runs up and down in III.D.4. 

Three examples of sample paths for different b values 
are given in Figures III. D. 5. 2 through III. D. 5. 4. The degen- 
eracy of the model is clearly present in the sample paths as 
a tendency to produce runs of equal, increasing or decreasing 
values, respectively. A comparison of Figures III. D. 5. 3 and 
III. D. 5. 4 quickly demonstrates that while these two sample 
paths have the same correlation, they produce significantly 
different (X n ) sequences. This is a graphic indication that 
non-normal processes are not determined solely by their 
correlation structure. 

Figures III. D. 5. 5 through III. D. 5. 7 are the scatter 
plots associated with the sample paths in Figures III. D. 5. 2 
through III. D. 5. 4, respectively. Here, too, the degeneracy 
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of the model is clearly present in the tendency of the 
(X n ,X n+1 ) plots to lie on the line X n+1 = bX n . The slope 
of this line determines whether the runs are of equal, in- 
creasing or decreasing value. 

6. Conditional Expectation and P (X .. > X ) for 
Antithetic Variables n n 

Results similar to those obtained in III.D.5 can be 
obtained for the moving minimum model with negative correlation. 
The procedure for determining the conditional expectations and 
the probability that X n+ ^ is greater than X n using antithetic 
variables is exactly the same as that in the previous section. 
First, the probability of each of the four possible combina- 
tions of X n and X n+ ^ values is computed, the conditional 
expectation or probability is computed for each case, and 
the final weighted sum of conditional expectations or proba- 
bilities is finally computed. In one instance no closed form 
answer is available and numerical procedures are used. 

Recall that the generation scheme when using antithetic 
variables is 

X = MIN([^]E, [b+l]E‘ ] , (III. D. 6.1) 

n d n n— l 

where (E n , n = 0,1,...} is an iid sequence of Exponentially 

distributed random variables with unit mean, (E', n = 0,1,...} 

n 

is generated from the (E } sequence by the relationship 
-E n 

E^ = - In (1-e n ) which implies that {E n } is also iid Exponen- 
tial with unit mean, b > 0. 
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First, consider the case where X = 



n 



[^~— ] E and 
d n 



X 



n+1 



f— — — 1 E 
L b J n+1 



Then 



P(X n = ^ E n' X n + l“ 



- p([ Tr 1E ni [b+1IE A-i' [b+11E A> 



Using the standard conditioning argument produces 



P (X„ = [— 1— ] E , X 



n 



n 



n+1 



= £±ii E 

L u J ^ 



n+1' 



00 

= / [b+l]E^, [^]E n+1 < -[b+l]ln(l-e" Y ) | E r 



y ) e Y dy 



= / P(E 

0 



n-1 



> 



£) P (E n+1 <_ -bln [1-e y ])e y ay 



/ e ~ y/b (1- [l-e~ y ] b ) e” y dy 
0 



r - (1+ b )y 

J e dy 

0 



00 k -w d+r) 

/ (1-e y ) b (e ' ) dy 



The first integral is straightforward. In the second integral, 
the change of variable u = e~ y , -- ^ = dy makes this integral 
recognizable as the integral of a Beta random variable. Making 
this change of variable and making the appropriate changes in 
the limits of integration produces 
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P(X = ,X = [b+1 ] E 

n bn' n+1 n+1 



- / (l-u) b (u) b du 
0 



b+1 



r (l+b) r (l+g) 

r ( 2+b+i) 



p (X n = t^i]E n< X n+1 . [b + l]E n+1 > 



b 2 r (b)r (|) 



b+1 (b 2 +b+l) (b 2 +l)r (b+i) 



(III. D. 6. 2) 



!_ l"| 

In the second case, X = [£?-=■] E and X , , = [b+l]E'. Pro- 

n on n+i n 

ceeding as before 



P(X = 



n 



[ Tr )E n' X n + l = lb+ 1 )E A> 



p HT 1E ni [b+ 1 IE A-l' lb+11E Ai 'X 1 E n*l' 



/ Q ?([ ^F 1] y - [t>+l]E^_ 1 ,-[b+l]ln[l-e" y ] < [^] E n+] _ | E n = y) e~ y dy 



< r b±I- 



/ P (E^_ 1 > P (E n+1 > -bln [1-e y ])e y dy 



/ e- y/b (l-e' y ) b e- y dy 
0 



00 k -d+cr )y 

/ (l-e - ) e dy 

0 
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This is the same integral as the second integral in the first 
case. Thus 



P(X = 



n 



[ TT 1 E n'Vl- > b+ 1 l E A J 



b 2 r(b) r(|) 



(b 2 +b+l) (b 2 +l) T (b+i) 



(III .D.6 . 3) 



Next consider the case where X = (b+l)E’ . and 

n n-1 

X n+1 - 'TT’W Then 



P< x n= (S±i)E n+1 ) 



= P( [b+1 ] E ’ < [^]E , [^]E ^ < [b+1 ] E ’ 
n — b n 1 b n+1 — 1 n 



/ P( [b+1 ] E^ _< [^p]y, [^p]E n+1 i [b+1 ] In [1-e y ] | E n = y) e _y dy 



/ p ( E n 1^) P (E n+1 £ -bln [1-e y ])e y dy 



oo oo _(l+i.)y 00 00 ( 1+fr) 

/ e y dy - / e ey - / e y (-e y )^*dy + / (1-e y )^(e y ) dy 

0 0 0 0 



The first two integrals do not present a problem. Making the 
change of variable z = 1-e y , dz = e y dy makes the third inte- 
gral easy. The last integral is the same as the second integral 
in the first case. So 
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P(X n = IWll^-Vi-lT 18 .*! 1 



= 1 - 



b 2 r (b) r (g) 



b+i b+l (b 2 +b+l) (b 2 +l) r (b+i) 



P(X = 



n 



[b+l ] E ' , , X , , = E . . ) 

n-1' n+1 b n+1 



b 2 r (b) r (i) 



(b 2 +b+l) (b 2 +l) T (b+J-) 

b 



(III. D. 6. 4) 



Finally, consider the case when X = (b+l)E' , and 

n n-1 



X 



n+1 = IMIIe;. Then 



P(X n = tb + l]E-_ 1( X n+1 - [b+l IE') 



- p ' ‘ b+1 1 E A i E „ < tb+i i e; £ c^i E n+ i > 



/ P([b+1]E^< [i±i]y,-[b+l]ln[l-e" y ] < C^]E n+1 |E n = y)e" y dy 



/ p(E n-iib )P(E n+l -bln [1-e y ])e" y dy 



/ (l-e- y )V y dy - / ( l-e -y ) b (e _y ) ~ dy 



d + i) 



0 



0 



These integrals are the same as the third and fourth integrals 
in case three. Therefore, 
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P (X = 



n 



Cb+ 1 ]E A-l' X n+l= [b+ 1 ]E A> 



b 2 r(b)r(^) 



b+l 



(b 2 +b+l) (b 2 +l)r(b+i) 



(III .D. 6 .5) 



The conditional expectations given a specific case 
of the values of X n and X n+ ^ can be written by inspection. 
Hence, 



Case 



x „ = lT IE n'Vl 



1 b J n+l' 



X = (~i]E n ,X = [b+l ] E ' ; 
n bn n+l n 



EjX n+1 iX n ^O 

b+l 
b * 

- (b+l) In (1-e b+1 ) 



x n = ' b+ 1 ' E A-l' X n + l= 



b+l 



X = [b+l ] E ’ ,,X = [b+l ] E ’ ; 

n n-1 n+l n 



b+l 



Combining these results with equations III. D. 6. 2 through 
III. D. 6. 5 and letting 



G = 



b 2 r (b) r (i) 



(b 2 +b+l) (b 2 +l) F (b+£) 



we have 
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(III. D. 6. 6) 



E(X n+l^ X n = y) = I E(X^ +1 \x n = y,case i)P(case i) 



i=l 



n+1 1 n 



L 

( b+T' G) ( ^T } " G (b+1) In (1-e b+1 ) 



+ G(^gi) + (g^-- G) (b+1) 



E(X n+ll X n = y) = 



2 - G (b+1) (1+ln [1-e b+1 ]) 



(III. D. 6. 7) 



Similarly, we can derive the expression for E ( x n I x n+ ]_ = Y) 



Case 



m n h= n+1 ^> 



X = E ,X . = . ; 

n b n n+1 b n+1 



b+1 
b * 



X = [^~]E , X , = [b+1 ] E ' ; 
n bn n+1 n 



__X_ 

,b+l N i b+1 ^ 
(-g— ) In (1-e ) 



X = [b+1 ] E ' . ,X .. = [~^]E , -i ; 
n n-1 n+1 b n+1 



(b+1) 



X = [b+1] E ' . , X , = [b+1 ] E ' ; 
n n-1 n+1 n 



(b+1) 



Then using III. D. 6. 2 through III. D. 6. 5 and again letting 



G = 



b 2 r (b) r (i) 



(b 2 +b+l) (b 2 +l) r (b+g) 



E ( x n I x n +i = Y) = 1 E (X n |X n+1 ,case i)P(case i) 



i=l 



n 1 n+1 



= (^) ( S | r -G)-G(^ln(l-e b+1 ) + (b+1) G+ (b+1 ) (g^g-G) 
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E< x n |x n+1 = y) 



(III .D.6 . 8) 



= 2 - G (— ^— ) (1 + In [1-e b+1 ]) 



The probability that X^ +1 is greater than X n can be 
approached in the same fashion as the conditional expecta- 
tions. The second case will be reserved for individual 
attention. 



Case 



P(X_+1 > XJ 
n n— 



X n = 



1 b J n+1 



1 

2 



X n- [b+11E A-l' X n + l= ^ 1E n + l 



b+1 



X = [b+1 ] E ' . ,X , = [b+1 ] E ' 
n n-l n+i n 



1 

2 



The second case, X = [^t-^JE and X . . 

n bn n+1 

allow a closed form solution. We get 



[b+1 ] E ' , does not 
n 



P(X n + l >X n> 



P( [b+ 1]E A > [tt±]E n ) 



— 7 ? 



P (- [b+1 ] In [1-e Jn ] > t^]E n ) 



-E E 

P (-In [1-e n ]>-^) 



-E b -E 
P ( [1-e n ] <e n ) 
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- y o b ' y o 

Hence, if we can find y Q such that (1-e ) = e , then 



-y 



0 



P(X n+ ^ > X n ) = 1-e . The required solution can be found by 

numerical means for any given value of b. A computer program 
to determine y^ to an accuracy of 10 for a given value of 
b and to compute P(X n+ ^ > X fl ) was prepared. A graph of the 
results is presented as Figure III. D. 6.1. When using anti- 
thetic variables, the moving minimum model has a restricted 
range of possible values for P(X n+ ^ > X n ) . The maximum value 
of approximately 0.509 occurs at about 0.30. The minimum 
value of approximately 0.491 occurs at about 3.33. 

This small range of values for the P (X . . > X ) is 

n+i n 

shown in the relative indistinguishability among the sample 
paths displayed in Figures III. D. 6. 2 through III. D. 6. 4. Of 
more interest are the scatter plots presented in Figures III. D. 6. 5 
through III. D. 6. 7. In these plots the degeneracy of the moving 
minimum model is clearly displayed. When X n achieves a value 
of y based on E , then X , , is constrained to have a value less 

-x n n+1 

than -In (1-e n ) . In the case where equality is achieved, the 

second case in the discussion in this section, the point (X , 

-X _x n +l 11 

X n+ ^) lies on the curve e n + e ‘ =1. This constraint 

is clearly evident in the scatter plots. Thus the moving 

minimum model displays a degenerate behavior for negative 

correlations as well as positive. 
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E. THE BETA-EXPONENTIAL MODEL 



1 . Introduction 

A third method that can be used to generate correlated, 
marginally Exponentially distributed random variables is a 
special case of the Beta-Gamma model given in Lawrance and 
Lewis [Ref. 18]. This model generates an {X n } sequence using 
the relation 



X n " + B n (1 -*'* )E n-l 



n = 1,2, 



(III. E. 1.1) 



where {B n (q,l-q), n = 1,2,...} is an iid sequence of Beta 
random variables, (B (l-q,q) , n = 1,2,...} is an iid sequence 
of Beta random variables, (E n , n = 0,1,...} is an iid sequence 
of Exponential random variables with unit mean, (B n (q,l-q)}, 
(B n (l-q,q)}, and {E n } are mutually independent, and 0 < q < 1. 
The density for a Beta(m,n) variable is 



r (m+n) m-1 .n-1 

r (m) r (n) (1 x; 



0<x<l;m>0;n>0. 



(III. E. 1.2) 



In practice the Beta random variables can be generated from 

two Gamma distributed random variables using the relationship 

B(m,n) = ; t 7 — r-T 7 T 7 — r from Kotz and Johnson [Ref. 19], where G(K) 
G (m; +G (n; 

is a Gamma random variable with a shape parameter of K and 
a scale parameter of one. 

This is a special case of the Gamma model considered 
in Chapter II of this thesis. It works because, as described 
by Lewis [Ref. 10], in III. E. 1.1 the product of the B^(q,l-q) 
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variable and the variable is a Gamma (q) variable. Simi- 
larly, the product of the B n (l-q,q) and the E n _j_ variables is 
a Gamma (1-q) , independent of the Gamma (q) variable. Conse- 
quently, their sum is an Exponential variable, X n . The {X n ) 
process is clearly one-dependent, as for the NEMA(l) process. 

Because of a lack of a closed form solution for the 
integral of the Beta density when the limits of integration 
are not from zero to one, this model is the least tractable 
of those considered in Chapter III of this thesis. However, 
its correlation structure can be determined, an expression 
for the Laplace-Stielt jes transform of a sum of r random varia- 
bles can be derived, and the probability of X n+1 being greater 
than X n can be simulated. An advantage of this model is that 
it extends directly to moving average Gamma processes (see 
Chapter II) . This extension is not possible with the NEMA(l) 
or the moving minimum model . 

2 . Correlation Structure, Positive and Negative 

The correlation structure can be determined using a 
standard approach. We have using III. E. 1.1 

Wl - (B n (q, 1 -q)V B n ( 1 -q,q)E n .i ) (B n+ i(q, 1 -q)E n+1 

+ B n + l ( 1 - q ' q)E n ) 

* B n+1 (q ' 1 - q> B n (q ' 1 - q) E n+l E n +B n+l (<3 ' 1_q) B n d-q,q) E^E^ 

+ B n+ 1 ( i-q , q) B n (q , i-q) E n +B n +i ( i-q . q> B n ( i-q . q) E n E n -i 
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Taking expectations and using the iid nature and independence 
of the (B n (q,l-q)}, {B n (l-q,q)}, and {E n } sequences yields 



E(X n X n+l } = g 2 + qd-q) + 2q(l-q) + (1-q) 2 



Hence, 



COV(X n ,X n+1 ) = q ( 1-q) 



and 



CORR(X n ,X n+1 ) = q(l-q), 0 < q < 1. (III. E. 2.1) 

As with the other linear additive models, this correlation is 
double valued and lies in the range (0,^-) . 

The range of possible correlations can be extended 
to negative values by modifying the generation formula by 
including an {E^} sequence. Thus 

X = B (q , 1-q) E_ + B (l-q,q)E’ ., (III. E. 2. 2) 

n n ^ n n n-1' 

where all random variables and constants are as defined for 
III. E. 1.1 and (E^, n = 0,1,...} is an iid sequence having a 
specified correlation with the {E n } sequence. In particular, 

E n and E^ may be an antithetic pair. The correlation of the 
{ X n } using II.E.2.2 can be determined in the same way as before. 
Consequently, 
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X X ^ 
n n+1 



(B n (q ' 1-q)E n + B n (1 ~ q,q) E n-1 ) 



(B n+1 (q,l-q)E n+1 + B n+1 (l 

= Vl'^-q) B n <qA-q> E n+1 E n + B 

+ B n+1 (l-q,q)B n (q,l-q)E n E^ + B 
Taking expectations as before 



-q#q) E ^) 

n+1 (q,l-q)B n 
n+1 (l-q,q) B n 



d-q,q)E n+ 1 E '. 1 

U-q.qlE'E^ 



E< x n x n+1 ) = q 2 +q(l-q)+q(l-q) [COV(E n ,E^)+lJ+(l-q) 2 

= 1 + q(l-q) COV(E n ,E^) 



Therefore, 



COV(X n ,X n+1 ) = q(l-q)COV(E n ,E^) 



and 



CORR(X n ,X n+1 ) = q (1-q) CORR(E n ,E^) , 0 < q < 1. (III. E. 2. 3) 

When E' = E , the correlation is one and III. E. 2. 3 
n n 

reduces to III. E. 2. 2. When E and E' are antithetic the 

n n 

correlation is -0.6449 and negative correlations result. When 

q is 0.50, the correlation for the {X R } sequence falls in the 

(-0.16,0.25) range depending on the correlation between E n 

and E ' . 
n 
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3 . Laplace-Stielt j es Transform of a Sum 



T = l X , (III. 3.1) 

r i=l 1 



where {X^} are defined by III. E. 1.1. Then 



T = X. + X 0 + . . . + X 
r 1 2 r 



= B 1 (q, 1-q) E 1 +B 1 (l-q,q)E Q +B 2 (q,l-q) E 2 +B 2 (l-q,q) E 1 



+ ... + B r (q,l-q)E r +B r (l-q,q)E r _ 1 



= B r (q, l-q)E r +B 1 ( 1-q , q) E Q + [Bj^ ( q , 1-q ) +B 2 (l-q,q) ]E 1 



+ • • • [ B r_l (q/ 1-q) + B r (1-q/q) ]E r _ 1 



-ST 

Let <}> T = E(e r ) . Then using the iid nature and independence 

r 

of {B (q,l-q)}, (B (l-q,q}, and { E } 

a 1 XT 1 1 



-S [B r (q, 1-q) E r +B 1 (l-q,q) E q + ^ B i (<3' 1-q) + B 2 (l-q,q) }E 1 

+ ... + {B . (q, 1-q) +B v .(l-q,q) }E ,] 

<1> T (s) = E ( e r_1 " r_ ' L 

r 

-sB (q, 1-q) E -sB (l-q,q)E 

= E(e r n )E(e 1 U ) 



-s [ ( B , (q,l-q)+B (l-q,q]E. 
x E(e 1 Z i ) 



“S (B . (q,l-q)+B (l-q,q) ]E . 

x . . .E(e r "' L r r X ) 
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t' t BE (s)) 



r-1 



r 




where 




E (e 



-S [B i (q,l-q)+B i _ 1 (l-q,q) ]E. _ i 



) . 



The Laplace transform of the sum of two Beta random variables 
is a confluent hypergeometric function. Its form is too compli- 
cated to be of significant value in deriving the analytic 
behavior of the Beta-Exponential model. 

4 . Empirical P(X n[ ^ >X n ) 



the probability of X n+ ^ being greater can not be analytically 
determined with a reasonable amount of effort. In an attempt 
to establish a range for this probability, a simulation was 
used. In order to achieve a precision of at least 0.001, 
sixty-eight thousand comparisons were generated for each of 
ten random number seeds. The Beta random variables were 
generated using the Kotz and Johnson [Ref. 32] relation 

s-* / __ \ 

B(m,n) = g~ (^j +g( n ) explained in III.E.l. The Exponential 
sequences were generated by a call to a standard generator of 
Exponentials. When the simulation was run for nineteen values 
of q from 0.05 to 0.9 5 in steps of 0.05, the P(X n+ ^ > ^ n ) was 
0.500 for all values of q. 



Because of the presence of the Beta random variables 
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Although the empirical probability that X n+ ^ is 
greater than X n is constant at a value of 0.500, reminiscent 
of the autoregressive model of Chapter II. B. 6, the distribu- 
tion of X n+ ^ -X is not symmetric and no simple proof for 
this situation has been found. 

The low serial correlation and the apparent invaria- 
bility of the p ( x n+ ]_ >x n ) makes the use of sample paths and 
scatter plots of little value. Samples are provided in Figures 
III. E. 4.1 through III. E. 4. 12. 
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FIGURE III. 
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FIGURE III. E. 4. 4 
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FIGURE III.E.4.6 
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FIGURE III. 
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IV. PRELIMINARY DATA ANALYSIS 



A. INTRODUCTION 

During the period 1955 through 1969 a weather ship sta- 
tioned in the Gulf of Alaska (50°N,145°W) collected, among 
other data, wind speed data at three hour intervals [Ref. 33]. 
The existence of the wind speed data was brought to Professor 
Lewis 1 attention when a student in Oceanography asked him to 
provide a model suitable for simulating wind velocity data. 

The simulated data was required as input to models of ocean 
temperature mixing. A copy of fifteen years of wind speed 
data was obtained for this thesis. The intent was to do a 
preliminary data analysis and then determine whether any of 
the models discussed in this thesis could provide an adequate 
representation of this data and, hence, a method for gener- 
ating wind velocity sample paths for oceanography simulation. 

Initially, the models discussed here are strong a priori 
candidates for data of this nature. Intuitively, there is a 
strong feeling that an assumption of independence among the 
data is unwarranted. Hence, autoregressive and moving aver- 
age models which are designed to reflect the behavior of 
correlated data should be considered likely candidates. The 
non-negative nature of the data mitigates against the use of 
classical time series anlaysis which is based on the assump- 
tion of a normal distribution, and hence negative values, for 
the series. The existence of zeros in the data tends to make 
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the use of transformations, like taking the log, somewhat 
less appealing than otherwise. These considerations indi- 
cate that the models discussed in this thesis should be 
considered as likely candidates for modeling the wind speed 
data. 



B. ANALYSIS OF THE RAW DATA 

There were 43,800 data points available for analysis, 

2920 for each of the fifteen years between 1955 and 1969 
inclusive (the extra data for leap years was discarded) . 

Since this size data base made it inconvenient, if not 
impossible, to manipulate by hand, each year's data was 
plotted as a means to promote familiarity with the data. 

The plot of each year's data and the plot of the data averaged 
over all fifteen years (e.g., all data taken at 0300 on 1 
January of each year were averaged) are presented in Figures 
IV.B.la through IV.B.lp. Several characteristics can be 
observed from these figures. First, and perhaps most obvious, 
is the expected yearly cycle of the data. Values at the 
beginning and end of the year tend to be higher than those 
in the middle. Second, the data is discretized to a large 
extent. There exist obvious horizontal lines of equal valued 
data. A sort and plot of the entire data set reveals that 
the data consists of values that are integral multiples of 
1.03 with a few values between these multiples. Next, on 
some occasions a series of high values will all be equal. 
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1956 Raw data 
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Ficjure IV. D. 1 c . 1957 Haw data 



m : 1 IUNU ! NG: NON FI 




233 



Piqure IV.B.Id. 1958 FUw data. 
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figure IV.B.le. 1909 Raw data 
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Figure IV.B.ly. 1961 Haw data. 
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1962 Raw data 
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1963 Raw data 
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1966 Haw data 
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Figure XV.B.lm. 1967 Raw data 
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Figuie IV.B.ln. 1968 Haw data 
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Figure IV.B.lo. 1969 Haw data 
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i.jmu IV. II. 1 |». J5 year average raw data. This sliows « he yoaily cyiAc < 1 < Y 



indicating that some clipping may have occurred (see Figure 
IV.B.la vicinity of 2575 and 2625, Figure IV.B.lb vicinity of 
400 and 525, etc.) . These last two characteristics indicate 
that statistical properties which are sensitive to the be- 
havior of the "tail" of a distribution may be affected. The 
final observation about the data is that there are apparently 
intervals when the data was not actually collected. These 
instances appear as reasonably long strings of values which 
have a strong linear appearance (as though the values were 
produced by linearly interpolating between two boundary 
values). See Figures IV.B.lh (vicinity 2400), IV.B.lj 
(vicinity 2250), IV.B.lk (vicinity 50 and 1750), and IV.B.lm 
(vicinity 150 and 2725) . 

The cyclical nature of the data is somewhat more apparent 
in the plot of the data averaged over the fifteen years (see 
Figure IV.B.lp). Additional evidence of this yearly cycle 
is presented in Figure IV. B. 2. This figure presents twelve 
box plots, one for each month. The data values plotted are 
the monthly average wind speed for each of the fifteen years. 
The interquartile range and extreme values are shown in a 
standard fashion. As an adjunct to this analysis of the 
year cycle, the coefficient of variation for the monthly 
averages was computed. The coefficient of variation was 
essentially constant. See Table IV.B.l. This will have 
an impact on the choice of the type of model used to model 
this data. 
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Table IV.B.l. Average wind velocity by month for each year. 
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This yearly cycle is also shown in the periodogram and the 

log of the periodogram of the data (averaged over 15 years) as 

presented in Figures IV. B. 3 and IV. B. 4, respectively. The 

periodogram is computed from the data in the following way. 

N 

Let {X , n = 1,2,...,N} be the raw data and let X = J X. be 

2 i=1 1 

the mean of the (X } sequence and a v be the variance. Let the 

n a 

(Y , n = 1,2,...,N} be formed from the { X } sequence using the 
n n * 3 

relation 



Y n = X n - X, (IV.B.l) 

where N = 2920 is an even number. The Fourier transform of the 

{ Y^} sequence will have both a real and complex component and 

will have j elements. Let (Z^, n = 1,2,...,^} be the Fourier 

transform of the (Y } sequence and let Z. D and Z. T be the real 

n jk 3 1 

th 

and imaginary components of the j element of {Z n }, respectively, 
t h 

Let Pj be the j element of the periodogram. Then 

Pj = (Z? R + Z^ x )/2-nNo^ f j = 1,2, ...,| (IV. B. 2) 

defines the periodogram of the { X n } sequence. 

The periodogram dramatically presents the yearly cycle 
(j = 1) as the dominant effect (P^ > 150) , although there is 
some indication of a six month cycle (j = 2, P 2 = 9.0). Some- 
what surprising is the apparent lack of any strong time of day 
effect. The log periodogram reinforces the dominant role of 
the yearly cycle and indicates that six month and six and 
twelve hour cycles (j =2, j = 1460, j = 2920 respectively) 
may be important. 
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The loy of the per i odo<| tarn of Fiqure IV. R. 3. shows the harmonies at 6 horns and 12 hours, and a backytuinui 
spectrum closely resemhl in«j that of a first-order aut oro«j ress i ve piocoss . 





The correlation structure of the data is presented in Table 
IV. B. 2. The column indicated as "15 Yr Avg" is the average of 
the values of the fifteen years. The column indicated as "15 
Yr SD" provides the standard deviation of the values about 
their mean. It is not the standard deviation of the average. 
This latter quantity can be obtained by dividing by the square 
root of 15. The last column provides the correlation structure 
of the average data. The estimated correlations remain artifi- 
cially high in this case because averaging reduces the varia- 
bility of the data about the year cycle which intensifies the 
artificial increase in correlation due to the year cycle. The 
correlation structure revealed in Table IV. B. 2 for individual 
years closely resembles that of an AR(1) model, in that the 

k-step correlation is approximately the one-step correlation 
th 

raised to the k power. The correlations in the table have 
a tendency to be slightly higher than the theoretical, calcu- 
lated value, but the agreement is reasonably good for about 
ten steps. Beyond that point the correlations are kept up by 
the year cycle, which is not as prominent in the yearly data 
as it is in the averaged data. If nothing else the disparity 
between the two correlations is evidence of the existence of 
a trend in the data. 

At this point sufficient information is available to de- 
termine some characteristics of the general form of the model 
for representing the wind speed data. As noted above, the 
correlation structure is similar to that of an AR(1) process 
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with positive correlation although nuances may appear as the 
year cycle is removed. Hence, an AR(1) model should be used 
as a starting point for the construction of the model for 
wind speed. 

In addition, the cyclic nature Of the process must be 
modeled. This can be done with either an additive or multi- 
plicative model. An additive model might have a structure as 
follows . 



n 



= y + e , 
H n n' 



(IV. B .3) 



where {X n } is the time series under consideration, y n is a 

deterministic function of n, and the innovative sequence 

{e n } is a stationary sequence of random variables. In the 

usual model this stationarity implies that the marginal vari- 
2 

ance a is constant and the correlations only depend on the 
lag (i.e., p( x n ' x n+ ] c ) = p (k) ) . Using the same definitions 
the multiplicative model would have the form 



X 



n 



y e , 

K n n' 



(IV. B. 4) 



where again the (e n ) sequence is stationary and independent 

of y . A characteristic of the additive model is that the 
n 

coefficient of variation is a function of the value of y . 

The multiplicative model produces a constant coefficient of 
variation. Since the data has a coefficient of variation that 
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is essentially constant, in the crude monthly analysis, the 
multiplicative model is preferred. 

We have yet to determine the exact form of the mean, 
in equation IV. B. 4. However, we do know that this mean will 
have a yearly cyclic nature. We also have yet to determine 
the general structural nature of the innovative process e . 
These subjects are addressed in the following sections. 

C. THE FORM OF THE MEAN; DETRENDING THE DATA 

Two basic models were considered to represent the mean. 

The first was a single harmonic sinusoidal model 

y n = a + b l sin ^2920^ + b 2 COS ^2920^ = a + k cos ^2920 + ' 

(IV.C.l) 

where k = (b^ + b^)^^ and 6 = tan b (- . The second was 

an exponential sine with one harmonic 



a+b. sin ( 



n 



= e 



27rn 

2920 



) +b-cos ( 



2irn 

2920 



) a kcos( :Sz7 +e) 

= e e 



(IV. C. 2) 



The second model has the theoretical advantage that it can 
not be negative and will represent higher harmonics in a com- 
pact form. The sinusoidal model may or may not be negative 
depending on the values for a, b^, and b 2 • In spite of the 
theoretical preference for the exponential sine, both models 
were used initially to see if either produced significantly 
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better results. Note that if k is small the models are 
hardly distinguishable. If k is large, the exponential sin 
will clip at low values and is a cycle that would have many 
harmonics in its Fourier transform. 

The values for the constants in equation IV.C.l were 
determined in a straightforward procedure using the least- 
squares regression procedure of MINITAB and the data averaged 
over 15 years. These estimates could also have been obtained 
from the periodogram at = 2 ttN, I(oj 1 ) ~ {(b 1 ) z + (b 2 ) z }, 
using the relations 

N x. 

a = l = X; (IV. C. 3) 

i=l N 

N 27ri 

b^ = 2 £ X^ sin ( - ■ ■) /N = imaginary component (IV. C. 4) 

i=l of periodogram at 

2tt/N ; 

N 27ri 

b 2 = 2 £ X . cos (— ^=-) /N = real component of 

i=l 1 periodogram at 2 tt/N. 

2 

The variance of these estimates is 2a /N if the X. 's are 

e i 

independent, but since this is clearly not the case here, 
estimates of the variance of the estimates cannot be obtained 
directly. The results of the estimation are contained in 
column 1 of Table IV.C.l. 

Similar results were obtained for the constants in IV. C. 2 
by a slightly more complicated procedure. In order to use a 
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TABLE IV.C.l 



Parameter Estimates for Models of the Mean Value 
Function of the Wind Velocity Data 

ESTIMATE 

1 harmonic 1 harmonic 2 harmonic 2 harmonic harmonic 

PARAMETER sine exp sine sine exp sine exp sine 



a/a' 10.230 2.309 10.230 2.307 



2.307 



b x -0.176 -0.011 -0.175 -0.011 -0.011 



b 2 2.560 0.260 2.566 0.260 0.260 



-0.593 -0.057 -0.057 



-0.397 -0.054 -0.054 



0.014 



0.001 



- 0.010 
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least squares approach, a linear relationship must be estab- 
lished for the mean value of the process. Taking logs is 
the obvious technique to employ, but this introduces a compli- 
cation. Taking logs and expectation of IV. C. 2 we have 

E(ln X n ) = In p n + E(ln e n ) 

, , . , 27m . , , , 27m . , 

a + b^ Sln ^ 29 20 + ^2 COS ^2920 + C ’ 

For example, if the {£ n l sequence is marginally distributed 
as a unit Gamma variate, G(l,k), then c = ^ (k) , where (k) 
is the digamma function (derivative of In r (k) ) . See Cox and 
Lewis [Ref. 29], pages 24-27. The value of the constant c 
will be combined with the constant a in the least squares 
estimation using the In X r 1 s , giving the constant a 1 = a+c. 

To estimate a+c without making Gamma assumptions for the 
innovative process, the X n 1 s are divided by 

^n b l sin 2920 + b 2 COS ^2920^ 

to give X^. The data is then divided by the average of the 

__ / 3 , ” 4*0 ) 

X ' s which estimate e . The result of this is a series 

n 

with mean value (within statistical fluctuation) of 1 if the 
model for the cycle is correct. The values obtained are listed 
in column 2 of TAble IV.C.l. The results of these estimates 
are in Figure IV.C.l. In this figure the average data 
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The disparity between the l wo-yoar trend functions and the data indicates the need for a second (6 month) 



is plotted against the value computed for u n using both of 
the models under consideration. 

When using a multiplicative model, the residuals are 
formed by dividing the raw data by the mean. The results of 
this procedure are presented in Figures IV.C.2a through IV.C.2p 
using the exponential sine model for the mean. The results 
are not significantly different using the sinusoidal model 
of the means. Hence, only the results for the average data 
are presented for this case in Figure IV. C. 3. 

The log periodogram of the average data detrended using 
the sinusoidal model for the mean is shown in Figure IV. C. 4. 

A five-step moving average of this log periodogram is pre- 
sented in Figure IV. C. 5. The detrending has clearly reduced 
the importance of the yearly cycle, but still shows some evi- 
dence of a six month cycle and six and twelve hour cycles. 
Similar information is provided for the average data detrended 
using the exponential sine model for the mean in Figures IV. C. 6 
and IV. C. 7. This model does not reduce the effect of the 
yearly cycle as much as the sinusoidal model for the mean, 

* but still shows the six month cycle as being important and 
some evidence of six and twelve hour cycles . 

Since the exponential sine has the theoretical advantage 
of being non-negative and both models of the mean produce 
similar results when applied to the data, the exponential 
sine is selected as the model of choice and the analysis is 
continued using it exclusively. 
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Figure* 1 V . C . 2l> . 1956 Detrended data 




1)11 I f) : Ml M I I ) l J H I 3- MCUJMI T MINI) \ 

ni. i nr no i NO: mul i i pi i cn 1 1 vi- i: xpqnt n i i m 




263 



Figure IV.C.2c. 1957 Detrended data 
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P'icjure IV.C.2d. 1958 Detrended data 
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Do trended data 
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Figure lV.C.2f. I960 Detrended data 
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Figure I V . C . 2 g . 1961 Detrended data 
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Figure IV.C.2U. 1962 Detrended data 
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Figure IV.C.2i. 1063 Detrended data 




J 
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Figure lV.C.2j. 1964 Detrended data 
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figure IV.C.2k. 1965 Do trended data 
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Fiqure IV.C.21. 1966 Detrended data 
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Figure lV.C.2m. 1967 Detrended da La 
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Fir,„ro lV.C.2n. 1968 fletromlp-l data 
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net rended data 
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Figure IV. C, 3. Average path detrended (sine) 
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Figure iv. C\ 4. hog per iodoqram shows dominance of six-month cycle and presence of six and twelve hour cycles after 

removal of yearly cycle. Background spectrum resembles first -order nut ore«ji ess i ve process. 






279 



i Htum Nt i 




i gg pi n i guggiu 

mini: HcsmiJfiL 3-nmim r winij vi kicitii 
IJI I m N l J I N G : Mill I I PL l CUT I VL L XPONT N 1 1 111 




280 



t icjuro 1V.C.6. I ,o<} per i odotjram shows dominance of six-month cycle and presence of six and twelve hour cycles aftei 

removal of yearly cycle. 
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Figure IV 










D. RESIDUAL PROCESS PROBABILITY STRUCTURE 

Having removed the dominant seasonal effect from the 
data, it is possible to investigate the structure of the 
residual process e n in the model 



X 



n 



V=n‘ 



( IV . D . 1 ) 



The two facets of the probability structure of the stationary 
process { e n } which were addressed in Chapter II were the 
marginal distribution and the correlation structure. The 
residuals produced by dividing the raw data by the appropriate 
value of the mean were supplied to HISTF, a histogram and box 
plot routine developed at the Naval Postgraduate School. 
Histograms for each year and the entire data set were pro- 
duced. These histograms are presented in Figures IV.D.la 
through IV.D.lp. The shape of the histograms is consistent 
over the years and indicates that a Gamma distribution is 
appropriate for modeling the innovative factors. The param- 
eter k can be estimated as the reciprocal of the coefficient 
of variation squared (see equation II. B. 4. 16) . The estimated 
value of k for each year is given in Table IV.D.l. 

A careful examination of the statistics associated with 
the histogram will reveal that the values for the skewness 
and kurtosis are low compared to the theoretical values for 
the Gamma distribution, namely 2//k and 6/k, respectively. 
However, this is not unexpected when one recalls the 



282 



FREQUENCIES SAMPLE SIZE • 2920 



o i 
i 
i 



l 

* 

i 

I 

I 

* 
^ i 
i 
* 



< r i 
<N I 



T> I 
n i 
I 

♦ 
* I 
I 

-» l 




I rn 
I 



* * 

* ♦ *J 

» * oo 
♦ — * 
♦ psi 



♦ * <M 

♦ * 

* — • 

♦ * <\l 

♦ * 

♦ 

* * * * 

* * * * B 
444* fO 

* * ♦ N 



4444 * 
44444 4" 

4*444 ^ 

44444444 

**♦*-4-*** 

44444444 

* 

4*44444 44 O 
444444444 J* 

4I»***444< -* 

4444444 444 
*-*♦*♦*♦*** 

4 

(444444444444444444 
(4444*444 4 44 44 44444 f- 



n | 
fN I 
I 



44444444444 
*44 444*4444 
***♦*♦**■*♦* 



4444444 
♦♦*4**4 
* 44 41 44 44 44 



— • I 44 ♦ 44 *** 44 ♦ 44 ♦ 44 44 44 * 44 ♦ * 44 44 44 * 44 44 44 ♦ 44 4 44 

r» | * 44 44 44 * * ♦ * 44 44 44 44 * * 44 ♦ * ♦ * « 44 44 44 44 * 44 ♦ * fO 

-* I 44 ♦ 44 44 44 * 44 44 -4 44 * 44 * 44 44 44 44 44 * * * 44 44 44 -4 44 44 * 



00 | 
00 | 
«“• I 



* 44 44 44 44 44 44-4 44 * 44 * 44 44 44* ***« **44*4*4-4*44 
**44«**4*****4*4****44*4*444*«44 



*44 44*4*44t«4*4 44444»«44*44 44*4 O' 

******* ♦♦♦♦**♦*♦♦****** ***** 4*4 -• 



N| ********* ♦ ***** * ** ************ 4 ****** * * -* 

<f 1 *4*44«4*«4 4444«44tl*4*«**4« 44**44444*44 

<M I ♦***♦♦***♦**********♦♦♦**♦***«********* 



CO I I I I I I I I I I | M I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 

oi • ••*•* 

mi *♦**♦♦♦*♦♦**»♦♦♦*♦******♦**♦♦**♦♦**♦***♦♦*♦***♦*** S' 

* ♦— • 
OI ***** ******4I ♦♦♦**♦♦* • ******************41* O 

m | ************ ***<*•»**.****•*-«**«*«**«**'»*'» 

<M| ****♦***♦»****-**************«♦********♦** 



O' I 

00 I 
I 



******* *********** ** *** *♦****♦***♦*♦♦*****»« 






********* *******4* ******** *♦* ******** 4 ♦♦*♦ 



******* *4** ♦-** **** ********* **** t> 
** I ******* *♦♦**#♦*« •** « < 4* *4* 4 * vf 



<*> I 

Nf I 

— I 



o 




u 



oooooo 

oooooo 

oooooo 

OOQfflO^ 

OJOUT*m*-aO 

o • • • 



* 

» 



-*z-» 

LU<UJ 

0—3 

— UJ— 

ISI 



* 

» 

* 

* 



5 

O' 

c 

TJ 

a 



* U 

Z. UJUJUJUJUJ * 4-> 

O JJJJJ * '!> 



.?Z2Z2 

r 

OOOOODD 

joooooi: 



oo ZomO'nox 

Q X ... * *x 



l/l — 1/4HH 

(— o -ONO 

Z » O' -I* I 
•jj uj nr»uj 

s » •O'oo 

O •OO'OO'T'O 

a: m>r -/(M 

—• O' OT' 

_j inr» in^ 

<t CDfVI 00 (\| 

Ct * • • • 



uj 

UJ l/>l/> 

i/l — 

X uji/> 

UJ zo^N 

X 2H<< 

O Jco-t- 

— » m^xoujuj 



in 

fM 

O 

in 

m f\i .oooo 

-•OONOO 

— • -o ~<oo 

OfM 0—0 
r»<M — <ooO 
<Nin vf »/*- 

. • * m . 



o 

uj x > < 

OXUJ UJ 
Q ZUJ>Q X 
< < O UJCL 

UJ — U.ZOV) 
X iOuj*ZQ 
a. «J(-auj< — 

1/1 >V0UJX£XX 



ooo—o 

ooono 
ooo-o in 
> •jO'-rro 
tj C'O'Otrur* 
o .O' . 0 s • 

Q 

UJ 



IU 

-I ZZO 

« Z4<Z 

a. -d LUUJ<| 

I- /-ITtt 

z *io—qo 

uj luiiia — — 

o i x»-rx 



3 

O 

JS 



I ^ 

| O' 

I *“• ' 



O -u 

I e 

I jj n> 

| o c 
l^o 
a a, 

X X 

◦ Cl 

-a 

o 

I T 5 > 
C H 
-0 JJ 
(0 

s o 
I <a -h 
w m 
31 a 

O *H 

U U 



o 

> 

M 

u 

3 



283 



FREQUENCIES SAMPLE SIZE - 2920 



•vi I 

I 



♦ TO 

* • 

* «>* 



♦ * A 
*— • • 
* * ♦ rvj 



® I * * * * * 

-* I *##♦♦ (*> 

• ♦**♦* H 

4 ♦ * 

r- I ******* <■* 

<vji +♦**+♦* 

♦ * 

0 I ******** 

fVJ I ******** A 

I + **»**♦♦ CT* 

♦ • 

f- I ********** H 

cnl 

I ********** 

* * 

01 ************* 

A I ************* r<- 

I ************* N 

♦ *— • 
N I ->*♦♦*♦*******♦**♦♦**♦♦*-< 

Q* I *********************** 

I **** *****♦♦**♦»**•***♦ * 

* * 

^1 *********************** 

3*1 *♦*♦*♦♦♦**♦♦♦**♦♦*♦•*♦♦* O 

I *********************** -3 

♦ • 

A | »* 4 *************************** * H 

I'll ******************************* 

— < I ********* *4****** *4 *4^4***^**** 

* ♦ 

fOI *«-t4**«4*4-»4*4*4**4****4**4l 

^ I ♦******♦****♦****♦******♦♦*♦ N 

-* i **************************** r 

mi ******************************** 

mi ******************************** 

—* i • 

OI *♦*****♦♦*'*♦*♦44*** ♦*♦****♦♦**♦♦**** *♦***♦♦* 

03 1 *** 4 *************************************** * *- 

I ******************************************** N 

* ♦ • 

— | ********************************************** r~i 

| ************* ************************* ******** 

-*l ***** ******************************* M ********* 

* * 

0 | ********************* ************************* 

<7*1 *•********>*******♦******•********♦*********** o 

-M ***** ******************* 4 **♦*********♦*♦♦♦*,».* * O 

* *— • • 

1^1 *********** *♦♦♦*♦*♦**♦***♦»♦****♦****♦****** * -4 

I I I I I I I I I « I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I | I I I I I I | I 

I • 

♦ * 

A | ********************************************* 

03 | * *******4 ♦ **♦♦*****■*-••***-■« ******************** A 

♦♦**'*♦♦*♦******♦*****♦ 4 ******** 4 *********** ♦ 03 

* *— • 

A | *********************** 4 **************** 4 ******** * O 

01 **************************** ********************** 

(Ml ************ ***************** ********************* 

* * 

A I ********************************* 

ml **«****-»**«»***********-»* **♦**♦*♦ -* 

I ********************************* N 

* *— • 

N I ******* 4 *♦*•«♦*♦»♦***♦***♦*******♦*♦******♦*♦**♦** ♦ o 

♦ ••••• 

I **************************** ********************** 

+ * 

•O I ************************************** ***** 

r- I ******************************************* m 

I *4 * ****************** * ************ v ***** 4 * * A 

* *— • 

Oil ♦*** ♦♦* ***** *♦** 4 * * * * * * 4 * * * * ♦ o 

I *******<***********«♦ + *«**,*♦ 

<-*1 ***************************** 

* * 

O | ********************************** 

mi *********************** •**+«•***** o 

— • | ***«**»«<*** ***********4**4* ****** m 

* ♦>— • 

33 I ******* I *****4* ************** o 

— | ******** **4**»****«k*»*«**»** 

| ***************************** 

A | ************** 

A | ************** 33 

I ************** •— 

* ■*> — • 

A | ********* O 

mi ********* 

I ********* 

CM . ♦♦*«#♦•» ***** 

sfl ************* 

I *•»***•**«*«** o 

r~ o a <■ m cm — • o 

O O O O O O O 



OOOOOO 

oooooo 

OOOOOO 
000*'fl* 
COAO'f’l-O'Ii 
om*>A • • * 



D 2ZZ ZZ 

n 

-* 0000030 

o' loaiooor 

A z. O AOAO X 

a r • • ♦ • 



A — 0~»— 

k- O -*030 

Z I flo I 

uj jj m — uj 

3T AO • •vl'P- 

O OmOO'f-* 

T. <T | oo 

O A OA 

_j — r>- — r- 



AA 

A^* 

UJ A 

ZC 3 -*N 



UJ <£>-►- 
m-r ^ oujuj 
acXA^oo 



TJ 

c 

<n 

S 

flj 



«-«00 
DNCNOO 
- 0-0 *MOO 
33 a> a-to 
003** 
ru -r ~r »r- 

• « .fVJ . 



Q 

> 



uj at; > <t 
OXUJ UJ 
ZUJ>Q n£ 
<Q UJ CL 

— • u. 2:0 a 

a: quj <zo 

< ►— OU)<- 

>viurai 



O' 

-H 

Cu 



-ooo— o 

VOOIMO 

NOAAO 

> 

o o cr cr' -r 

4* . sAO'Of' • 



UJ I 

_J 7 ZO | 

< Z< <Z I 

UL <UJUJ< * 

)- z — ttx 

Z <0—00 
UJ UJUi<X— — 
o ZInXX 



284 



♦♦♦00000 0000000004 



FREQUENCIES SAMPLE SIZE * 2920 



o 



O I 
I 
I 

♦ 

°l 

I 

♦ 
o I 
I 
I 

+ 

o | 

I 

l 

<n i 



I 

I 

I 



oooooo 

oooooo 

oooooo 

OOONf'N 

mo 

o^-ow* • • • 



* * 
#* 
» * 



UJ<UJ 

0—0 



O I 
O I 
I 

♦ 
o I 

o I 
I 



♦**♦*4 

**♦**+*♦ 

*4*»4*4« 

*♦♦»♦*+* 

***♦♦■#*♦♦** 

4*4444*4444 

4*4*444444* 

M44*4 4*4 1 4444 
»444*4444*444* 



* * 
* * 
* * 



* ♦ 

+ - 
* ♦ 



P» I 
O I 
I 

4 

>f I 

00 I 
I 

4» 

— 1 

2 I 

♦ 

mi 44* 444 

I 4**44* 

— ' I * H4 ♦♦♦♦ 4 ♦+♦**■**♦♦* H* -#■*+♦ -♦ ■**■** ■* m 

■4 ♦— • 

'♦ I ***4444**4*44**44444#*4*44 44 4#44#4444*44* H 

— I 44444444444 4 44*4444 4 444»44*444*44*4444444 

fM I 

"* ♦ 

01 444**4#*#*44*#44*#**44*44#4*#44*4* 

N | 44444*4444444444**4444444*4*444444 m 

"• I 444*444444444**4 *44*44 44444*444444 — 

* + — • 

'O I *4444*44444444444444*44444444444444*444 -* 

Ol 444 4*4 44 4 44 44444*44 4««444*444444444*444 

<N I 444444444 4*4 44 4****44 *4*4444*4*4 44*44*4 

■* 4 

0 I 44*44444444444««4444*4 44*44*4444444444444444* -O 

^ I 4r ♦ 4 4 *** 4 4 *#♦♦* H* ♦+•*♦♦♦**++♦**♦***■* 4 4- ♦*.* * O 

+ • 1 1 « 1 1 1 1 • » « 1 1 1 1 1 t 1 1 1 1 1 1 1 1 » 1 1 1 1 r 1 • 1 r 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 — • 

N I O 

I ♦♦♦♦♦♦4 1444*444444*44444444 ♦♦♦♦♦♦♦♦♦♦ 

* ■*• 

m| 444 44 44 4444 4444*44 ) 444 44 44 44 4444 444 4*4 4 *444 4 44 

sT I *44444«»44»*4*44*«#4444 4444444«»*4«4 444*44*444 

(M I 444«4«444*t«4««4*44444444*4444444*«*4444444444 f*» 

*■ *■— * 

01 444444444444444«444I44I444444«4^*4 C 

r*l 444« 44i*44444«4*44l4 44«4 *4<«44->4«44 

— • I *4it4ir«44444*444444444444444444«44 

•*> ♦ 

mi 4444444444« 4 444 4 «4 4 4 144 444444444 444*44444*44 44444* 

>01 © 

CU | 444444 4 4 444 4 4*4 *4***444 4 444 4444 4 44 4 4 444 4 4 44 4 444 444 m 

T I 4444»44«44«I 44 4444444444*44444444*I4»4I4>4 O 

'M | *4444411*44 44*««44«44«(4«*«44«44 44lt444 4*» 

Nl *•***■*♦.«.*♦♦•«*♦«'•**-«*■**•*.«**•**’•■•*♦*■•♦♦••»«** 

♦ * 

m| 444*44*4*444*44*44|4444444|4«4 

mi 44 44*4444*44*^44444444444444 »» O 

— • I 4-4 4444444*44*4444444444l44t444.© 

4 4— • 

<N| 44444444*444*44444 O 

0»| ♦ *♦♦*■*•«♦■•■*•**■«■**•** 

| «44444<4444I44I44* 






O 22^22 

ca 

— 13^333 

<x z joojor 

4— — —a 

m 20momox 

Q 



m — 

>-» o o»mo 

2 I — I 
UJ j J’O'JJ 
2! 0>0 • • — >0 
O vfOOOm — 
z om i rsj -r 
>orj >Of\i 

-J Om o m 



mm 
m— 
aim 
.rO—rsl 
24- << 
-U0C4— 4- 
m>f*:0UJUJ 

Z srm^aa a3 



mr*» •mo- 
rn oomoo 
o>r moo 

— »CT* —1^0 

— m r-o-o 
rsivf m »>o 

• • »m • 



UJ QC> < 
0 ><lUJ UJ 
O ^UJ>0 cC 

< < C x a. 

UJ — < iZUH 
a£ afQuj«JZQ 
a. <h-cuj«j — 
m >i^uza: X 



mOO— O 
NOO<UO 
crccmm 

— C-UNfl 

>o mmmrn 

m v crm • 



0 
r— 4 

a. 

* 

o 

xi 



444444* 

• ♦♦*44* 

• 4*4^4* 



7<<JZ 

<iu UJ «1 

7-rr« 

•to—oo 

UJ UJO.*— — 

TP-ir 



285 



***+ *0000000000000+ 



FREQUENCIES SAMPLE SIZE - 2920 



°*I 

l 



O I 
I 

+ 
>#• I 



m l 
I 
I 

♦ 
o I 
I 
I 

♦ 
<Nl I 
-* I 
I 

♦ 

(M I 
<M I 

I 



******* 



<M I ********* 

'I’l **♦*♦***♦ 

I **♦***■#♦* 

♦ * 

mi ********** 

-t i ***♦*♦*•#** 

i * * * ******* 

♦ *- 

-* i * -*♦*■»****♦-»■*♦ * 

*• 1 *************** 

i «*******»*«*** 

+• * 

rg | ****** •****««*»** 

03 | *»**»*****■»*«***« 

I * + *■**■*•*•**■****■**♦* 

♦ +- 

rO I ♦♦*******♦♦**♦**♦♦***♦ 

-* I ******♦*******+*•****** 

-* I ***•**********♦**♦*♦*♦* 

* 

m i ****************** * *** * 

•*•*•*♦**♦«******#** + *•*** 

— • i *****♦*♦*♦♦ *********** 

♦ +- 

oi *♦****♦♦**♦**+ ********-*********** 

(*• | »**************^******»«********* 

^ ^ ******#■#***# + •*#*#****♦•*******•**** 

mi ♦ *•#*♦*-*♦* + * + *#******* *****+■♦**♦* 

OI A******************************* 

-•I ******,***********+♦************ 

♦ + - 

r i + *♦**■***-**•***♦**•**•***** + ♦ + ***■* 

m I #**«-r*** ************** ******** 

— I *♦******♦**♦***♦♦♦♦*♦******+*+ 

(Ml ******** **♦♦*♦*♦***♦* *************** *********** 

<r \ ********************** >»♦** *♦*♦♦*** ***♦♦*♦*♦*♦** 

<M I **J********************** ******** ************** 

* *— 

mi I I I | I I | I i i I i i I i I i I i i I I I I i I i i i i | I I I | i i i i i i i i i 

0 i ************************************ 

I 

* * 

!> I #*-»*♦**♦**♦♦******♦«**** ** 4***********4 *********** 
mi *********** ***»»***«*«•• * ******** *»■**♦* ************ 
<M I ******************* ******** ******** ♦♦*<****♦**♦*** 

* +*” 

CD 1 ************************************* 

031 *****«■**♦■* 4 ** 44 *** « * * * * *♦*«*♦*♦♦***♦* 

-• I **<******^«»**4 ♦******«♦♦*♦♦*«*****♦* 

♦ * 

01 * * * »**************************«*************** 

mi ********************************************** 

(Ml *#**************•«•♦** **•«**♦* ******************* 

+• 

i ********* **♦■**♦♦**♦** ********* ******** *■#***♦***♦ 

■t i ***♦*♦»***♦♦** ******** *************** * ********** 

(Ml ****■**************»#**♦♦«*****»♦*** + ******♦****♦ 

* * 

OI ***4* ******* «»♦<******♦«*♦***** 

mi ****»+*****««*»«* ************** 

i **********♦*«*******<**«****«** 

* * — 

O I *************************** 

mi **************** *********** 

-* i *************************** 

* * 

vT I ********************* 

OI ****** •«»***•♦***♦**** 

— I ********************* 

* +■ — 

mi ************* 

oi ***♦*«*♦♦*«** 

i ************* 

* ♦ 

— i *********** 

mi *********** 

i •*•««****** 



a> 



m 



m 



(M 



oooooo 
oooooo 
oooooo 
ooovroN 
r'*— mr*-m 
omom • • • 
♦ • • •-*-*(<> 



UI<U» 

0—0 

2QZ 



Z. UIUJUIUIUJ 



•o 

a 

CC 



ZZZZZ 

ooooooo 

?OOOOOS 



m zomOmox 

— * — -»i\imf>*o , < 

a 



u 

u 

o 

'O 



m — mo-* 
i— O omo 
Z I *-0 | 
ai ui UVMUI 
X — *(M • •vrJ* 

o moom vf 
X (Mm m- -» 
mm ■rm 

~i -*vf — *vT 

<x oo«m a>(M 




Z 

UI 

o mm 

to- 
es luo 

ui zq-in 

X 

o uicSl— t— 

— m^^jouiui 

X 



<r 

(M 

m 

— »o •-too 
<tt MOOOO 
mo — « o o 
m(M o'imo 
( s»(M -*mm 
r^m ~r «r>- 



a 

x 

o 

J3 

C 

KJ 

S 

•-3 

» M 



T3 

-H 

Q 

> 

M 



o 

ui ae> < 
o ><uj ui 
a ziu>o ct 

< <t O UlCL 
ui — LLZOO 
ac aCQUKZQ 
a <i— Ouj< — 
m xoozaJ 



v* 

2 

O' 

Cm 



f^OO(MO 

ooooc 
ojOmmc 
> oor-mo 
vj O'Tmmo 

Z •U'O'O' • 

UI -* • • • — 

o 



I- 



o 



UJ 

zzo 

OUUJ *X 

7-rza: 

-1 Q-QU 

u>uj*x— — • 
ni-zi 



286 



++ ♦+ +++ * + +** 0000000000000 + 



FREQUENCIES SAMPLE SIZE - 2920 



l 



® I 
I 
I 

O I 
— I 
I 

«\ I 

I 



<NI I 

<N | 



O I 
cm I 
I 



* * 

***•♦ + <M 
***** 

***** <M 
***** 
***** 



I 

* I 
I 



*4 * ♦* *•** INI 
******** 
******** 



****** 



► **•♦**■*** CO 

♦ — • 

► ■**■*.**♦♦* — 



>A | 
-t I 
-< I 



► ***-♦•*-'1’*** <M 



^*4»*****»* 



^ I 

H | 
<NI I 
♦ 
<r i 

CO I 



I 

CO I 

CM | 



m I 

rO I 
<M I 



♦ ** ** — 

* * * * * -* 



ft******************************** 
•*•**♦•#*>*■*<•*♦*♦♦*♦*♦♦**•#■•♦*-* •*■♦**•»**■*■**♦*••** 
I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 1 I I I I I I I I I I 



**>* * *** * *•* * ******* **** * *** **** ******** *** O 

► **■*■#** •*♦♦** **-*•*•* * * * * * * •***■»-• *•#♦*-***■***■***■** 
»*«*•*«***** ** *« ******* * **« ************ * ** 

M * *♦* * ****** * * ♦****♦**♦♦* * **l « *•*•***>**« ** 
****** ********************l************** O 
*•■*****■**•*■*•«**•*■**■**«#♦•**■*** •*•*■«***♦****■♦<•#** P- 



***»***«*«< 



h* 4 * *41 



********i 



**•*.-***>** * * 



► •*•*****• 



m « 
— I 
-• l 



41**444*144* 

* 

******* 



OOOOOO 

oooooo 

OOOOOO 

OOOOOH 

o>r -o • • • • 

• • • —4-4-44*1 



0—0 

204 



UJUJUJUJUJ 



o zzzzz 

o 

- 3D33333 

o C I03000Z 

co zOinO’riOx 

Q 



0»f- 

-•O 

J»iA 

r-O* 

— **o • • a? m 

O O' 00^0 
<M<M ~KP 

O' <7" CP CQ 

r»>C0 OcO 

-4fA — «o 



U1<S) 

UJi/i 

20-*(\l 

ujac»-*- 

o sT JUJUJ 

3 ; 3i^^.aco 



X> A *iAOO 

cr-ocnoo 

'J'fM OOO 



oj x> < 

U>'tUJ UJ 

a z uj >o <x 
< <tO aiQ. 
uj — u. 2 TO*/> 

04 C£QUJ<2Q 

a. • 

v/> 



inooco 

OCN^m 
P~0— — ^A 
OOOOco 



zro 
Z <1 z 

vlULU< 

r-rrtt 
•4C— oo 
ujil'oc 1 — 

ZJ'-TT 



7 



1 



287 



Figure IV.D.le. Histogram and boxplot of 1959 detrended data 
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FREQUENCIES SAMPLE SIZE - 2920 
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Figure iV.D.lo. Histogram and boxplot of 1969 detrended data 
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Table IV. D . 1 



Moment Estimate of the Gamma Shape Parameter by Year 



Year 


Estimate 


1955 


3.96 


1956 


4.16 


1957 


4.38 


1958 


3.68 


1959 


3.65 


1960 


4.45 


1961 


3.87 


1962 


3.87 


1963 


3.86 


1964 


4.49 


1965 


4.32 


1966 


4.04 


1967 


5.08 


1968 


4.24 


1969 


5.41 
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discretization and clipping that has apparently occurred 
while the data was collected and processed. As a rough 
check on the extent of the clipping required to produce 
values for the kurtosis and skewness similar to those for the 
data the following procedure was examined. A sample of 2920 
values for a Gamma distribution were produced with mean 1 and 
shape parameter k = 4 . 0 by a call to the NPS random number 
generator LLRANDOMII (SGAI1A) . A histogram of these values was 
produced sequentially for the following cases. The data was 
clipped so that all values over four were set equal to four, 
all values over three were set equal to three, and finally the 
highest ten percent of the data was set equal to the value of 
the 2890 sample order statistic. The first four central moments 
were estimated under each of the conditions . The results are 
presented in Table IV. D. 2. The results indicate that a 
clipping of the top ten percent of the data will yield results 
for skewness and kurtosis comparable with those observed in 
the data. 

TABLE IV.D.2 





Mean 


SD 


CV 


Skewness 


Kurtosis 


Gamma 


0.986 


0.504 


0.511 


1.128 


2.114 


Cut at 4.0 


0.986 


0.504 


0.511 


1.128 


2.114 


Cut at 3.0 


0.984 


0.498 


0.506 


0.994 


1.165 


10% Cut 


0.982 


0.489 


0.498 


0.861 


0.516 
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